dynamic shear experiment

ORIGINAL RESEARCH article

An experimental study on scale effect in dynamic shear properties of high-damping rubber bearings.

• Bridgestone Corporation, Seismic Isolation and Vibration Control Products Development Department, Yokohama, Japan

High-damping rubber bearing (HDRB) is one of the most popular devices used for seismic isolation of structures. In order to clarify the mechanical characteristics of HDRB, various loading tests have been conducted on the bearings and obtained data have been applied to practical design of isolation systems. In this study, in order to investigate the scale effect on the physical characteristics of HDRB, dynamic loading tests were conducted with full scale and scaled model isolators, which have diameters of 1,000 and 225 mm, respectively. The test program covers shear strain dependence tests, frequency dependence tests, and repeated loading dependence tests. Special attention is paid in the differences of shear characteristics caused by the specimen scale. Repeated loading test was conducted only with a scaled model, and the relationship of temperature increase of the specimen and shear characteristics was evaluated. In parallel, finite element analysis (FEA) of the isolator under repeated loading was conducted. After the constitutive model of FEA was identified by the results, the FEA was extrapolated to simulate repeated loading of 1,000- and 1,600-mm-diameter isolators, which cannot be tested realistically by dynamic loading. Change of properties along the increasing number of cycles and temperature distribution of full scale and scaled down isolators were investigated. Necessity of consideration for the scale effect in the evaluation of HDRB properties by dynamic testing is discussed.

Introduction

Seismic isolation technology has gained popularity in the recent decades as one of the measures for seismic protection of structures ( Murota, 2009 ; Nishi and Murota, 2012 ). Seismic isolation is an aseismic design concept to reduce the seismic force transmitted to the structure by supporting it with a flexible support member to elongate the natural period of the structure and thereby decouples it from the ground effects. Basically, seismic isolation systems provide functions of restoring force and energy dissipation. The seismic isolation bearing (SIB), made up with layers of alternating rubber and steel plates, is the most popular device for providing restoring force and damping characteristics. In addition to conventional SIB, innovation of new types of SIB has been progressed by many researchers ( Losanno et al., 2019 ; Madera Sierra et al., 2019 ).

Many kinds of SIB testing, such as shear strain dependence, temperature dependence, frequency dependence, vertical loading force dependence ( Tubaldi et al., 2016 ; Kalfas et al., 2017 ), aging effect ( Hamaguchi et al., 2009 ), and ultimate properties ( Nishi et al., 2019 ), have been conducted and useful data have been obtained for the design of the isolation systems. In order to evaluate realistic isolator characteristics, full scale isolator testing under dynamic loading is desired ( Infanti et al., 2004 ; Yamamoto et al., 2012 ; Kato et al., 2014 ). However, in the case of dynamic loading of such specimens, a substantial capacity of hydraulic systems including the number of accumulators, high-level control devices, and high-precision measurement systems are required for the test setup. Therefore, full scale testing is generally conducted under quasi-static loading conditions, and rate-dependent characteristics are evaluated using a scaled model by a relatively small capacity dynamic testing machine. As an example, in ASCE 7–16, chapter 17 ( American Society of Civil Engineers [ASCE], 2017 ), it is permitted to use a scaled model for the test of isolators that has rate-dependent properties, considering practical reasons caused by the limitation of the number of facilities, which provides sufficient capacity for dynamic tests of full scale specimens. In ISO 22762 “Elastomeric seismic-protection isolators” ( ISO 22762, 2018a , b ), the allowable scaling for each test type is specified. Basically, rubber material itself has frequency dependency in its restoring force characteristics. High-damping rubber bearing (HDRB) is one of the types of SIB that has relatively large frequency dependence in terms of shear properties. Generally, the rubber material of HDRB is filled with carbon or other ingredients, and interaction occurring between the polymer and the filler under stretching condition dissipates kinetic energy as heat builds up. Increasing temperature of rubber results in the decrease in shear stiffness and damping ratio. The thermal conductivity of rubber is much lower than that of metals, and the build-up heat in rubber radiates from the surface of SIB. Therefore, the accumulation of heat inside rubber is affected by the scale of length, which is considered as the main reason of scale effect on the physical properties of SIB.

In this study, firstly, the scale effect on the fundamental shear properties of HDRB is investigated by dynamic loading tests with two types of specimen: full scale and scaled-down specimens with a scale ratio of 1/4.44. Specimen diameters are 1,000 and 225 mm. The measured properties are shear strain dependence and frequency dependence. Secondly, dynamic repeated loading test with 200 cycles is conducted on isolator specimens with a diameter of 225 mm.

After the 2003 Tokachi-Oki earthquake, long-period and long duration ground motions have become one of the most highlighted topics in earthquake engineering in Japan. A long duration ground motion may continue for over 5 min. and in the case of a seismically isolated structure, isolation devices are subjected to vibration for a long period ( Building Research Institute Japan, 2016 ). On this background, this study deals with the subject of scale effect on isolator characteristics under repeated cyclic loading. During the repeated loading, the properties of isolators are changed according to the increasing number of cycles. This property change firstly occurred by loading history ( Tubaldi et al., 2017 ) for first several cycles, and then the increase in temperature is mostly influenced by the large number of cycles. During the tests, the surface temperature of the specimens is measured, and the relationship between the change of mechanical properties and the temperature according to the increase in the number of cycles is investigated. Additionally, finite element analyses (FEAs) are carried out under identical conditions with the tests, where parameter identification is conducted. After confirming the accuracy of the FEA model, the dynamic repeated loading characteristics with 1,000- and 1,600-mm-diameter isolators are investigated with this model and compared with the test results.

Dynamic Loading Test

Definition of shear properties.

The shear characteristics, such as effective shear stiffness K eq , equivalent shear modulus G eq , equivalent damping ratio H eq , and dissipated energy E d , of HDRB are defined in Figure 1 . K eq is defined as the slope of the straight line from the point of (maximum load, maximum displacement) to the point of (minimum load, minimum displacement). Under sinusoidal loading, the force at maximum displacement has generally small difference from the maximum force. In the definition of the viscoelastic characteristics, the force at maximum displacement is generally used for the estimation of effective stiffness. However, the maximum displacement point, where the shear force drops sharply, is very difficult to identify and sometimes causes loss of accuracy in the calculation of characteristics when the time step for data acquisition is not precise enough. For this reason, effective stiffness is defined as mentioned above in this study. In frequency dependence test and shear strain dependence test, shear properties at the third cyclic loop are used for evaluation.

Figure 1. Definition of shear properties of HDR.

Test Specimens

The dimensional characteristics, performance specifications, and material properties of full scale and scaled-down models are shown in Table 1 and Figure 2 . The full scale and scaled model isolators have consistency of important parameters such as unit rubber layer thickness, number of layers, and first and second shape factors. However, because of the manufacturing process, there are slight differences for each parameter, which are indicated in Table 1 . For reinforcing plates and flanges, it is challenging to follow an exact scale. Furthermore, full scale isolator has an inner hole with a diameter of 25 mm for manufacturing purpose. The authors considered these differences in the scale as negligible in this study. Full scale and scaled model isolators are identified as scale-I and scale-II, respectively. Two different types of high-damping rubber materials are used in the specimens. They are identified as rubber-A and -B, which have equivalent shear modulus G eq of 0.392 and 0.620 MPa at 100% shear strain, respectively. The equivalent damping ratio H eq of both rubber materials is 24% at 100% shear strain. These two high-damping rubber materials had been improved in terms of load history dependence. A stress softening behavior, known as the Mullins effect, is improved compared with conventional high-damping rubber materials ( Murota et al., 2007 ).

Table 1. Dimensional characteristics and material properties of test specimens.

Figure 2. Test specimens.

Each test specimen is named according to size, rubber material, and specimen no. as follows: [size: I or II]-[rubber: A or B]-[No.: 1, 2, or 3].

Shear strain γ is defined as X max / h r , where X max is the maximum horizontal displacement at a cycle of shear loading and h r is the total rubber height in an isolator unit. The shear strain of 100% corresponds to the shear displacement equivalent to the total rubber height. Each test specimen was manufactured from a different lot of rubber material and vulcanization conditions. Therefore, the variation of the modulus of used rubber material and the cure state of rubber by vulcanization may affect the shear properties of isolators. Especially, the cure state of scale-I and scale-II may have significant difference in network structures of polymer, sulfur, and fillers. This may be also considered as a scale effect.

Test Facility

The dynamic loading tests of the scale-I specimen were carried out by “Seismic Response Modification Device (SRMD) Test Machine” ( Seible et al., 2000 ) in the University of California at San Diego. The loading plate of the test device, where the isolator is installed, slides over low-friction hydrostatic bearings. The generated friction force during dynamic loading was deducted from the raw data during data analysis. The friction force was measured by the UCSD laboratory in advance. The inertia force is calculated by acceleration and weight of mobile parts of the test device and was also deducted from the raw data. All test procedures were operated by the staff of the SRMD test facility in the UCSD.

The tests of scale-II were carried out by dynamic loading test machine in Bridgestone Corporation Technical Center, Yokohama, Japan. The maximum vertical load is 1,000 kN, and the maximum horizontal load and maximum strokes of the horizontal actuator is 200 kN and ± 0.3 m, respectively. The load cell, for measurement of shear and compression load, was installed just beneath the specimen. Therefore, the shear force obtained by the load cell does not include the friction force generated in the slide guide of the test machine. The inertia force is included but it is considered as negligibly small.

Test Conditions

Correspondence between the specimen number and the test type is shown in Table 2 along with test frequency and no. of cycles incorporated. Benchmark test was conducted to investigate the fundamental performance of isolators at a shear strain of 100% under a compressive stress of 13 MPa for material A and 15 MPa for material B, which is the test condition of the isolator used for determining the nominal shear modulus and the damping ratio by the manufacturer ( Table 2 ). In frequency dependence tests, the benchmark tests were conducted repeatedly between each test at specific frequency in order to evaluate the influence of loading history on the results, which should be properly considered at evaluation. The change of the properties measured by the benchmark tests is an indicator of the fatigue by numerous loadings to the specimen. Ambient temperature during the testing of the scaled model was controlled at 20 ± 5 degree Celsius. The temperature of scale-I in the UCSD was not controlled, but the measured temperature ranges between 20 and 26 degrees Celsius during the test. Therefore, the authors consider the influence of ambient temperature to the properties as insignificant and negligible.

Table 2. Summary of test conditions.

In frequency dependence test, the shear strain is 100%, the wave form is sinusoidal, and the number of cycles is three. A frequency of 0.33 Hz is considered as the standard vibration frequency in this study. The series of frequency is 0.01, 0.033, 0.1, and 0.3 Hz. Benchmark test was conducted between each test at a specific frequency as prescribed. In this study, frequency dependence is defined as the ratio of G eq , H eq , and E d at each frequency level to those of 0.33-Hz frequency.

In the shear strain dependence test, its series is 10, 25, 50, 100, 150, 200, and 270%. The test frequency is fixed as 0.33 Hz, and the wave form is sinusoidal. The benchmark test is conducted at the beginning and end of the test series. In this study, shear strain dependence is defined as the ratio of G eq , H eq , and E d at each shear strain to those of 100% shear strain.

In the repeated loading test, the specimen was subjected to a compressive load according to the corresponding nominal stress of 13 and 15 MPa for materials A and B, respectively. Under the controlled compressive load, the cyclic loading was conducted in the shear direction for a shear strain of 200%, with a frequency of 0.33 Hz. The total number of cycles was 200. The cumulative displacement was 72 m, which corresponds to 320 m for a full scale isolator unit assuming a total rubber height of 200 mm. The ratio of shear properties, G eq and H eq , at each cycle was computed by normalizing each value with that of the third cyclic loading. During testing, the surface temperature of the specimen was measured by a radiation thermometer.

Verification of Initial Properties of Test Specimens

Initial properties of test specimens were evaluated by the first set of benchmark tests. Comparing the nominal value of each property, all results were within the range of ± 20%, concluding that all test specimens were properly manufactured within designated margins. In the tests, scale-I specimens gave higher stiffness and a lower damping ratio than those of scale-II specimens. Figure 3 shows a comparison of the shear stress-strain relationship of I-A-1 and II-A-1, I-B-1, and I-B-1, for all three cycles. The factors that have possible influence on the difference in properties between scale-I and -II are the difference in the rubber material lot, vulcanization conditions, and the test device. Here, the difference in stiffness may be caused by the rubber material. The difference in the dissipated energy is quite small between scale-I and scale-II specimens, while the difference in stiffness is larger, which is attributed to the difference of the material lot and/or vulcanization conditions.

Figure 3. Shear stress-strain relationship of scale-I and-II of rubber-A and -B for all three cycles.

Although there is some difference in measured properties between scale-I and -II specimens, all specimens satisfy the standard deviation, ± 15% of the design value, which is generally considered as the acceptance criterion of isolators in practical use.

Results of Frequency Dependence and Shear Strain Dependence Tests

The comparison of shear stress–strain relationships under frequency of 0.01, 0.033, 0.1, and 0.33Hz for scale-I and -II with rubber-A and -B at the third cycle is shown in Figure 4 . Shear properties, G eq , H eq , and E d at each frequency, normalized by the value at 0.33 Hz are shown in Figure 5 . As it was expected, the shear modulus becomes higher for increasing values of the frequency. The shear modulus at 0.01-Hz frequency, which is considered as quasi-static, is more than 25% lower than that at the 0.33-Hz level. The damping ratio and dissipated energy values also show similar tendency. The difference in the frequency dependence between scale-I and -II is relatively small. Especially, for E d , results almost agree with each other. These results indicate that the test results of the scale-I isolator conducted by quasi-static test conditions can be corrected to the results of dynamic test conditions by applying frequency dependence obtained by dynamic test of the scaled model.

Figure 4. Shear stress-strain relationship of scale-I and -II in frequency dependence test at 3rd cycle.

Figure 5. Comparison of shear properties in frequency dependency of-Scale-I and -II at 3rd cycle.

Figure 6 shows the transition of each property in benchmark tests conducted before each test at each frequency. The first benchmark test [BT(1), number in () : set number of the test] is conducted for verification of the initial performance of the specimens as prescribed in 2.5. The total number of sets was 5, and their sequence is

Figure 6. Change of characteristics in benchmark test [BT(1) to BT(5)] of scale-I and -II at 3rd cvcle.

The absolute value of each BT was normalized by the value at BT(1). All properties in each test specimen decrease as the number of sets increases. Also, it can be observed that scale-II shows more reduction of stiffness and constant reduction in dissipated energy. However, correlation with scaling is not obvious. The change ratio varies approximately between −10 and −20%. The results indicate that during the prototype test of isolators, effects of fatigue by accumulated loading should be adequately considered when making judgments according to design criteria. It is suggested that when the frequency dependence test will be conducted, fatigue effect should be measured by the benchmark test between each test as conducted in this study, and adequate correction of the results should be carried out such as deduction of the decreased amount of properties by fatigue from the results.

The shear stress–strain relationships for all cycles and each property normalized by the value at a shear strain of 100% for both rubber-A and -B with scale-I and II at the third cycle of each shear strain are shown in Figures 7 , 8 . Significant difference in the absolute value of shear stress is observed in the shear stress–strain relationship for scale-I and II, which is shown in Figure 7 . The authors consider that the difference in stress is caused by the variation of rubber material and vulcanization process in scale-I and II specimens. However, it is noteworthy that the normalized value shows good agreement in scale-I and -II as indicated in Figure 8 . Deviation in shear modulus, ± 10 to 20% for example, does not affect the shear strain dependence. This result suggests that the shear strain dependence can be effectively evaluated using scaled models instead of full scale similar to frequency dependence.

Figure 7. Comparison of shear stress-strain relationship in shear strain dependence test of scale-I and -II for all three cycles.

Figure 8. Comparison of shear strain dependency of scale-I and -II at 3rd cycle.

The results in frequency and shear strain dependence tests show that even when there is significant difference in the absolute value of shear properties in scale-I and -II, the normalized trends with frequency and shear strain amplitude were not affected by the scale. In both tests, loading was conducted in three cycles. The temperature on the rubber surface on scale-I isolators increased only 4 to 5 degrees Celsius at the end of the test compared to the initial condition. Benchmark test results indicate effective stiffness and dissipated energy decrease as loading experience increases. Appropriate consideration should be made in the evaluation of the results in continuous sets of loading.

Results of Repeated Loading Test

Figure 9 shows the shear force–displacement relationship of each specimen in repeated loading tests for the entire cyclic loading protocol. It is observed that the shear force and E d are decreased by the increased number of cycles. The change of G eq , H eq , and E d under repeated loading is indicated in Figure 10 . The change of properties by repeated loading is considered as a combination of fatigue and temperature effect. For the first few cycles, fatigue effect, also called as “Mullins’s effect” ( Mullins, 1969 ) is dominant in the change of properties. After a few cycles, properties are majorly affected by temperature increase. At the 200th cycle, G eq decreased for approximately 40% in II-A-3 and 30% in II-B-3. There is no significant difference in the change ratio of G eq and E d for both rubber-A and -B. The increase in the surface temperature along the increase in the number of loading cycle, and the relationship between volumetric E d of rubber in specimens ( V r ) and surface temperature for both rubber materials are shown in Figure 11 . The temperature increases from 20 degrees to a maximum of 80 degrees Celsius in the case of specimen B, which has higher shear modulus and energy dissipation. The result indicates that the specific heat capacities of rubber-A and -B have no significant difference. Therefore, the difference of temperature increase in A and B is simply caused by the difference in accumulated energy dissipation during repeated loading.

Figure 9. Shear stress-strain relationship of scale- II-A and -II-B under repeated loading for shear strain 200% × 200 cycles.

Figure 10. Relationship of change of properties and number of cycles in scale-II-A and -II-B.

Figure 11. Change of temperature by increase of loading cycle and E d / V r in scale-II-A, and -II-B.

Finite Element Analysis for Evaluation of Scale Effect in Repeated Loading

As repeated loading tests on full scale specimens are quite challenging in existing test machines, FEA was conducted to predict the scale effect on the repeated loading by extrapolation to a large size isolator test case. FEA code named as the “Deformation History Integral Type (DHI)” model ( Mori et al., 2012 ; Masaki et al., 2017 ) was implemented for the study. The function of the code is a heat-mechanics coupled analysis, which consists of hyper elastic stress–strain analysis and heat-transfer analysis. The constitutive law of the mathematical model ( Mori et al., 2012 ) involves parameters for temperature affect and fatigue affect. Using both parameters, the properties of isolators during repeated loading are reproduced. The conceptual flow of the heat-mechanics coupled FEA and the comparison of the test results and FEA results for scaled model II-B-3 are shown in Figure 12 .

Figure 12. Concept of heat-mechanics coupled FEA and comparison of shear stress-strain relationship of Testing (Scale II-B) and FEA results.

Heat transfer boundary is set surrounded by isolators, and the coefficient of heat transfer was identified so as to agree with the surface temperature of the model. As explained in a previous publication ( Mori et al., 2012 ), parameters were identified by the results of scaled models and shear-block specimens. Firstly, the parameters for temperature effect on the change of properties were identified by the test conducted with shear-block specimen, which was conditioned in constant temperature under −10, 0, 10, 20, 30, and 40 degrees Celsius. The number of loading cycles was three, and the properties at the third cycle was measured and recorded. Property change in the temperature range over 40 degrees Celsius was assumed by extrapolation with curve fitting. Parameters related to fatigue were identified using data of property change obtained during repeated loading dependency tests with 200 cycles.

Using these parameters, the repeated loading for full scale isolators with a diameter of 1,000 mm, which was used in this study as specimen I-B, and 1,600-mm-diameter was simulated by an FEA model, and the results were compared with II-3. The analysis was conducted based on 50 cycles. The dimensional characteristics of 1,600-mm isolators are as follows:

• outer diameter = 1,600 mm, inner diameter = 80 mm

• thickness of unit rubber layer = 10.4 mm, number of lamination = 19, total rubber thickness = 198 mm

• thickness of reinforcing plate = 5.8 mm, rubber material = B.

Figure 13 shows a comparison of normalized analysis results for G eq , H eq , and E d with respect to the number of loading cycles, where values are normalized according to the results of the third cycle. The temperature distribution in all three isolator models at the final step is also shown in Figure 13 . The results show the scale effect on the dynamic characteristics under repeated loading, which was not observed in frequency and shear strain dependence tests conducted for the three cycles. Analyses have shown that the effect of increasing repeated cycles on the shear modulus is almost the same for each diameter. However, in terms of dissipated energy E d , there are significant differences between scaled model II-B and the largest size of the 1,600-mm-diameter isolator as the number of cycles increases. Obviously, the difference comes from the different temperature levels inside isolators, as the heat generation by high-damping rubber is proportional to the cube of the size (volume), whereas heat dissipation from the surface of the isolator is proportional to the square of the size (surface area).

Figure 13. Comparison of change in G eq , H eq , and E d of isolators with diameter of 225 mm (scale-II-B). 1000 mm (scale I-B), and 1600 mm in repeated loading by FEA.

The scale effect on the dynamic shear properties of HDRB was investigated by dynamic loading tests on full scale specimens with a diameter of 1,000 mm, and scaled model with 225 mm, using two different types of rubber material with soft and hard shear moduli.

Firstly, the nominal shear properties at 100% shear strain were measured and the manufacturing conformity of all specimens was verified. Then, dynamic loading tests of frequency and shear strain dependence were conducted. In both types of rubber materials, although there is some degree of difference between shear stress–strain relationship of full scale and scaled specimens, no significant difference was observed in neither dependences.

During the frequency dependence test, the change of fundamental properties was evaluated by the benchmark test, which was conducted between each test at a specific frequency. The stiffness and dissipated energy vales were decreased as the number of test cases increased. This fact suggests that when the continuous test is conducted, such as the prototype test of an actual project, the fatigue condition of the specimen should be appropriately considered in the evaluation of the results.

Next, repeated dynamic loading for 200 cycles of 200% shear strain with a scaled isolator specimen was conducted, and the relationship between the change of shear properties, temperature increase, and the number of loading cycles was investigated. As the number of cycles increases, the temperature of the isolator increased and stiffness and dissipated energy decreased.

The FEA model, which was developed in a previous study by authors for heat-mechanics coupled analysis, was implemented in order to investigate the cyclic characteristics of 1,000- and 1,600-mm-diameter isolators. The parameters in the constitutive law of the FEA was identified with test results of the scaled model with the diameter of 225 mm. The results show a significant difference in the change of shear properties by the number of repeated cycles. The results indicate that when the isolator is subjected to repeated loading over three cycles, the scale effect on shear properties is significant.

It is concluded that the scale effect on shear properties under a limited number of cycles, such as fewer than 10 cycles, can be neglected. It can be fed back to the practical case of the isolator test. In the prototype tests, if the frequency dependence is evaluated in advance, testing with a scaled model can be accepted. However, when long duration seismic input is considered, which has been a current issue since recent major earthquakes in Japan, the isolators may be subjected to an extreme repeated number of shear loading cycles. Considering these cases, the scale effect on the dynamic properties should be properly considered.

In this study, investigation of the scale effect is limited to the properties between 100 and 200% shear strain. The scale effect of HDRB isolators on ultimate properties such as shear breaking or buckling is under consideration as future subjects. Furthermore, investigation of the scale effect on other types of isolators, especially lead-core rubber bearing (LRB), is also considered in the next step. Energy dissipation is concentrated in the lead core where the thermal diffusion is considered to be affected by the dimension of the isolator.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest

NM and TM were employed by the company Bridgestone Corporation.

Acknowledgments

The authors would like to express great appreciation to all the members of the University of California at San Diego for their contribution to the dynamic loading test of full scale isolators.

American Society of Civil Engineers [ASCE], (2017). “Chapter 17 seismic design requirements for seismically isolated structures,” in Minimum Design Loads and Associated Criteria for Buildings and Other Structures American Society of Civil Engineering , Bombai: American Society of Civil Engineers.

Google Scholar

Building Research Institute Japan, (2016). “Study on repeated loading characteristics of seismic isolation devices and its influence to seismically isolated structures,” in Building Research Report, No.170 , 2016.4, Japanese.

Hamaguchi, H., Samejima, Y., and Kani, N. (2009). “A study of aging effect on rubber bearings after about twenty years in use,” in 11the World Conference on Seismic Isolation, Energy Dissipation and Active Vibration Control of Structures , Guangzou, 7–21.

Infanti, S., Papanikolas, P., Benzoni, G., and Castellano, M. G. (2004). “Rion-antirion bridge: design and full scale testing of the seismic protection devices,” in 13th World Conference on Earthquake Engineering , Vancouver.

Iso 22762, (2018a). “Elastomeric seismic-protection isolators,” in Part-1 Test methods International Organization for Standardization , Geneva: ISO.

Iso 22762, (2018b). “Elastomeric seismic-protection isolators,” in Part-3 Applications of buildings - specifications, International Organization for Standardization , Geneva: ISO.

Kalfas, K., Mitoulis, S., and Katakalos, K. (2017). Numerical study on the response of steel-laminated elastomeric bearings subjected to variable axial loads and development of local tensile stresses. Eng. Struct. Col. 134, 346–357. doi: 10.1016/j.engstruct.2016.12.015

CrossRef Full Text | Google Scholar

Kato, H., Akutsu, S., Mori, T., Nakamura, M., and Murota, N. (2014). “Development of large-scale elastomeric seismic isolation bearings for buildings,” in 5th Asia Conference on Earthquake Engineering , Taiwan.

Losanno, D., Dierra, I. E., Spizzuoco, M., Marulanda, J., and Thomson, P. (2019). Experimental assessment and analytical modeling of novel fiber-reinforced isolators in unbounded configuration. Composize Struct. 212, 66–82. doi: 10.1016/j.compstruct.2019.01.026

Madera Sierra, I. E., Losanno, D., Srano, S., Marulanda, J., and Thomson, P. (2019). Development and experimental behavior of HDR seismic isolators for low-rise residential buildings. Eng. Struct. 183, 894–906. doi: 10.1016/j.engstruct.2019.01.037

Masaki, N., Mori, T., Murota, N., and Kasai, K. (2017). “Validation of hysteresis model of deformation-integral type of high damping rubber bearings,” in 16WCEE, proceedings of 12th World Conference of Earthquake Engineering , Santiago.

Mori, T., Kato, H., and Murota, N. (2012). Modeling for analysis of high damping rubber bearings with effects of fatigue and temperature increase due to repeated deformation. Struct. Constr. Eng. 77, 1517–1526. doi: 10.3130/aijs.77.1517

Mullins, L. (1969). Softening of rubber by deformation. Rubber Chem. Technol. 42, 339–362. doi: 10.5254/1.3539210

Murota, N. (2009). Earthquake protection materials – reviews and future directions of elastomeric isolators. Soc. Polym. Sci. 58,

Murota, N., Nishimura, K., Kikuchi, T., Okutsu, N., Kuwamoto, K., and Suzuki, S. (2007). “Development of new high-damping rubber bearings with improved material for load-history dependence,” in Summaries of Technical Papers of Annual Meeting Architectural Institute of Japan , Japanese, 879–880.

Nishi, T., and Murota, N. (2012). Elastomeric seismic-protection isolators for buildings and bridges. Chin. J. Polym. Sci. 31, 50–57.

Nishi, T., Suzuki, S., Aoki, M., Sawada, T., and Fukuda, S. (2019). International investigation of shear displacement capacity of various elastomeric seismic-protection isolators for buildings. J. Rubber Res. 22, 33–41. doi: 10.1007/s42464-019-00006-x

Seible, F., Benzoni, G., Filiatrault, A., Post, T., and Mellon, D. (2000). “Performance validation of large seismic response modification devices,” in Proceedings of 12th World Conference of Earthquake Engineering , Auckland.

Tubaldi, E., Mitoulis, S., Ahmadi, H., and Muhr, A. (2016). A parametric study on the axial behaviour of elastomeric isolators in multi-span bridges subjected to horizontal excitation. Bull. Earthqu. Eng. 14, 1285–1310. doi: 10.1007/s10518-016-9876-9

Tubaldi, E., Ragni, L., Dall’Asta, A., Ahmadi, H., and Muhr, A. (2017). stress softening behavior of HDNR bearings: modelling and influence on the seismic response of isolated structures. Earthq. Eng. Struct. Dyn. 46, 2033–2054.

Yamamoto, M., Minewaki, S., Yoneda, H., and Higashino, M. (2012). Nonlinear behavior of high-damping rubber bearings under horizontal bidirectional loading: full scale tests and analytical modeling. Earthq. Eng. Struct. Dyn. 41, 1845–1860. doi: 10.1002/eqe.2161

Keywords : seismic isolation, high-damping rubber bearing, full scale, dynamic loading, scale effect

Citation: Murota N and Mori T (2020) An Experimental Study on Scale Effect in Dynamic Shear Properties of High-Damping Rubber Bearings. Front. Built Environ. 6:37. doi: 10.3389/fbuil.2020.00037

Received: 22 November 2019; Accepted: 11 March 2020; Published: 28 April 2020.

Reviewed by:

Copyright © 2020 Murota and Mori. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY) . The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Nobuo Murota, [email protected]

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Dynamic Experiments using Simultaneous Compression and Shear Loading

• Published: 10 July 2017
• Volume 57 , pages 1359–1369, ( 2017 )

Cite this article

• B. Claus   ORCID: orcid.org/0000-0002-0351-9866 1 ,
• M. Beason 2 ,
• H. Liao 1 ,
• B. Martin 3 &
• W. Chen 1  

717 Accesses

5 Citations

Explore all metrics

Material characterization at high strain rates under simultaneous compression and shear loading has been a challenge due to the differing normal and shear wave speeds. An experimental technique utilizing the compression Kolsky bar apparatus was developed to apply dynamic compression and shear loading on a specimen nearly simultaneously. Synchronization between the compression and shear loading was realized by generating the torsion wave near the specimen which minimizes the time difference between the arrival of the compression and torsion waves. This modified Kolsky bar makes it possible to characterize the dynamic response of a material to combined compression and shear impact loading. This method can also be applied to study dynamic friction behavior across an interface under controlled loading conditions. The feasibility of this method is demonstrated in the dynamic characterization of a simulant polymer bonded explosive material.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Subscribe and save.

• Get 10 units per month
• Download Article/Chapter or eBook
• 1 Unit = 1 Article or 1 Chapter
• Cancel anytime

Price includes VAT (Russian Federation)

Instant access to the full article PDF.

Rent this article via DeepDyve

Institutional subscriptions

Similar content being viewed by others

An Apparatus for Tensile Characterization of Materials within the Upper Intermediate Strain Rate Regime

A Novel Apparatus and Methodology to Characterise the High-Rate Behaviour of Materials Under Complex Loading Conditions

Inertial and Frictional Effects in Dynamic Compression Testing

Chen WW, Song B (2010) Split Hopkinson (Kolsky) bar: design, testing and applications. Springer Science & Business Media

Kolsky H (1949) An investigation of the mechanical properties of materials at very high rates of loading. Proc Phys Soc Sect B 62(11):676

Article   Google Scholar  

Christensen R, Swanson S, Brown W (1972) Split-hopkinson-bar tests on rock under confining pressure. Exp Mech 12(11): 508–513

Frew DJ, Akers SA, Chen W, Green ML (2010) Development of a dynamic triaxial kolsky bar. Meas Sci Technol 21(10): 105–704

Duffy J, Campbell J, Hawley R (1971) On the use of a torsional split hopkinson bar to study rate effects in 1100-0 aluminum. J Appl Mech 38(1):83–91

Nemat-Nasser S, Isaacs JB, Starrett JE (1991) Hopkinson techniques for dynamic recovery experiments. Proc R Soc Lond A Math Phys Sci 435(1894):371–391

Baker WE, Yew C (1966) Strain-rate effects in the propagation of torsional plastic waves. J Appl Mech 33 (4):917–923

Espinosa H, Patanella A, Fischer M (2000) A novel dynamic friction experiment using a modified kolsky bar apparatus. Exp Mech 40(2):138–153

Rajagopalan S, Prakash V (1999) A modified torsional kolsky bar for investigating dynamic friction. Exp Mech 39(4):295–303

Claus B, Nie X, Martin B, Chen W (2015) A side-impact torsion kolsky bar. Exp Mech, pp. 1–8

Nie X, Prabhu R, Chen W, Caruthers JM, Weerasooriya T (2011) A kolsky torsion bar technique for characterization of dynamic shear response of soft materials. Exp Mech 51(9):1527–1534

Ogawa K (1997) Impact friction test method by applying stress wave. Exp Mech 37(4):398–402

Article   MathSciNet   Google Scholar  

Huang H, Feng R (2004) A study of the dynamic tribological response of closed fracture surface pairs by kolsky-bar compression-shear experiment. Int J Solids Struct 41(11):2821–2835

Huang H, Feng R (2007) Dynamic friction of sic surfaces: a torsional kolsky bar tribometer study. Tribol Lett 27(3):329–338

Lewis JL, Goldsmith W (1973) A biaxial split hopkinson bar for simultaneous torsion and compression. Rev Sci Instrum 44(7):811–813

Rittel D, Lee S, Ravichandran G (2002) A shear-compression specimen for large strain testing. Exp Mech 42(1):58–64

Zhao J, Knauss W, Ravichandran G (2009) A new shear-compression-specimen for determining quasistatic and dynamic polymer properties. Exp Mech 49(3):427

Zhao P, Lu F, Chen R, Lin Y, Li J, Lu L, Sun G (2011) A technique for combined dynamic compression-shear test. Rev Sci Instrum 82(3):035–110

Lin Y, Qin J, Lu F, Chen R, Li X (2014) Dynamic friction coefficient of two plastics against aluminum under impact loading. Tribol Int 79:26–31

Download references

Acknowledgements

This research was supported by a cooperative agreement between Air Force Research Laboratory and Purdue University (FA8651-13-2-0005). We also acknowledge for the inspiration from Dr. Dan Casem of Army Research Laboratory in the design and application of this experimental technique.

Author information

Authors and affiliations.

Purdue University, 701 W. Stadium Ave, West Lafayette, IN, 47907, USA

B. Claus, J. Chu, H. Liao & W. Chen

Purdue University, 500 Allison Road, West Lafayette, IN, 47907, USA

U.S. Air Force Research Laboratory, 101 W. Eglin Blvd. Suite 135, Eglin AFB, FL, 32542, USA

You can also search for this author in PubMed   Google Scholar

Corresponding author

Correspondence to B. Claus .

Rights and permissions

Reprints and permissions

About this article

Claus, B., Chu, J., Beason, M. et al. Dynamic Experiments using Simultaneous Compression and Shear Loading. Exp Mech 57 , 1359–1369 (2017). https://doi.org/10.1007/s11340-017-0310-2

Download citation

Received : 14 March 2017

Accepted : 30 June 2017

Published : 10 July 2017

Issue Date : November 2017

DOI : https://doi.org/10.1007/s11340-017-0310-2

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

• Simultaneous torsion and compression
• Split Hopkinson pressure bar
• Find a journal
• Publish with us
• Track your research

Dynamic Shear Rheometer

The dynamic shear rheometer (DSR) (Figure 1 and Figure 2) is used to characterize the viscous and elastic behavior of asphalt binders at medium to high temperatures. This characterization is used in the Superpave PG asphalt binder specification. As with other Superpave binder tests, the actual temperatures anticipated in the area where the asphalt binder will be placed determine the test temperatures used.

The basic DSR test uses a thin asphalt binder sample (Figure 3) sandwiched between two circular plates. The lower plate is fixed while the upper plate oscillates back and forth across the sample at 10 rad/sec (1.59 Hz) to create a shearing action (Figure 4). DSR tests are conducted on unaged, RTFO aged and PAV aged asphalt binder samples. The test is largely software controlled.

The standard dynamic shear rheometer test is:

• AASHTO T 315: Determining the Rheological Properties of Asphalt Binder Using a Dynamic Shear Rheometer (DSR)

Asphalt binders are viscoelastic. This means they behave partly like an elastic solid (deformation due to loading is recoverable – it is able to return to its original shape after a load is removed) and partly like a viscous liquid (deformation due to loading is non-recoverable – it cannot return to its original shape after a load is removed). Having been used in the plastics industry for years, the DSR is capable of quantifying both elastic and viscous properties. This makes it well suited for characterizing asphalt binders in the in-service pavement temperature range.

The DSR measures a specimen’s complex shear modulus (G*) and phase angle (δ). The complex shear modulus (G*) can be considered the sample’s total resistance to deformation when repeatedly sheared, while the phase angle (δ), is the lag between the applied shear stress and the resulting shear strain (Figure 5). The larger the phase angle (δ), the more viscous the material. Phase angle (δ) limiting values are:

• Purely elastic material: δ = 0 degrees
• Purely viscous material: δ = 90 degrees

The specified DSR oscillation rate of 10 radians/second (1.59 Hz) is meant to simulate the shearing action corresponding to a traffic speed of about 55 mph (90 km/hr).

DSR Superpave Specification Logic

G* and δ are used as predictors of HMA rutting and fatigue cracking. Early in pavement life rutting is the main concern, while later in pavement life fatigue cracking becomes the major concern.

Rutting Prevention

In order to resist rutting, an asphalt binder should be stiff (it should not deform too much) and it should be elastic (it should be able to return to its original shape after load deformation). Therefore, the complex shear modulus elastic portion, G*/sinδ (Figure 6), should be large. When rutting is of greatest concern (during an HMA pavement’s early and mid-life), a minimum value for the elastic component of the complex shear modulus is specified. Intuitively, the higher the G* value, the stiffer the asphalt binder is (able to resist deformation), and the lower the δ value, the greater the elastic portion of G* is (able to recover its original shape after being deformed by a load).

Another way to look at this is that rutting is basically a cyclic loading phenomenon. With each traffic cycle, work is being done to deform the pavement surface. Part of this work is recovered by the elastic rebound of the pavement surface, while part is dissipated in the form of permanent deformation, heat, cracking and crack propagation. Therefore, in order to minimize rutting, the amount of work dissipated per loading cycle should be minimized. The work dissipated per loading cycle at a constant stress can be expressed as:

W c = work dissipated per load cycle

σ = stress applied during load cycle

G* = complex modulus

δ = phase angle

In order to minimize the work dissipated per loading cycle, the parameter G*/sinδ should be be maximized. Therefore, minimum values for G*/sinδ for the DSR tests conducted on unaged asphalt binder and RTFO aged asphalt binder are specified.

Fatigue Cracking Prevention

In order to resist fatigue cracking, an asphalt binder should be elastic (able to dissipate energy by rebounding and not cracking) but not too stiff (excessively stiff substances will crack rather than deform-then-rebound). Therefore, the complex shear modulus viscous portion, G*sinδ (Figure 5), should be a minimum. When fatigue cracking is of greatest concern (late in an HMA pavement’s life), a maximum value for the viscous component of the complex shear modulus is specified.

Another way to look at this is that fatigue cracking can be considered a stress-controlled phenomenon in thick HMA pavements and a strain-controlled phenomenon in thin HMA pavements. Since fatigue cracking is more prevalent in thin pavements, the parameter of most concern for fatigue resistance can be considered a strain-controlled one. With each traffic cycle, work is being done to deform the pavement surface. Part of this work is recovered by the elastic rebound of the pavement surface, while part is dissipated in the form of permanent deformation, heat, cracking and crack propagation. The lower the amount of energy dissipated per loading cycle the less likely fatigue cracking is. Therefore, in order to minimize fatigue cracking the amount of work dissipated per loading cycle should be minimized. The work dissipated per loading cycle at a constant strain can be expressed as:

ε 0 = strain during load cycle

This relationship between G*sinδ and fatigue cracking is more tenuous than the rutting relationship.

In order to minimize the work dissipated per loading cycle, the parameter G*sinδ should be minimized. Therefore, maximum values for G*sinδ for the DSR tests conducted on PAV aged asphalt binder are specified.

Test Description

The following description is a brief summary of the test. It is not a complete procedure and should not be used to perform the test. The complete test procedure can be found in:

A small sample of asphalt binder is sandwiched between two plates. The test temperature, specimen size and plate diameter depend upon the type of asphalt binder being tested. Unaged asphalt binder and RTFO residue are tested at the high temperature specification for a given performance grade (PG) binder using a specimen 0.04 inches (1 mm) thick and 1 inch (25 mm) in diameter. PAV residue is tested at lower temperatures, however these temperatures are significantly above the low temperature specification for a given PG binder. These lower temperatures make the specimen quite stiff, which results in small measured phase angles (δ). Therefore, a thicker sample (0.08 inches (2 mm)) with a smaller diameter (0.315 inches (8 mm)) is used so that a measurable phase angle (δ) can be determined. Figure 7 shows major DSR equipment.

Test temperatures greater than 115°F (46°C) use a sample 0.04 inches (1 mm) thick and 1 inch (25 mm) in diameter, while test temperatures between 39°F and 104°F (4°C and 40°C) use a sample 0.08 inches (2 mm) thick and 0.315 inches (8 mm) in diameter (Figure 8). The test specimen is kept at near constant temperature by heating and cooling a surrounding environmental chamber. The top plate oscillates at 10 rad/sec (1.59 Hz) in a sinusoidal waveform while the equipment measures the maximum applied stress, the resulting maximum strain, and the time lag between them. The software then automatically calculates the complex modulus (G*) and phase angle (δ). Much of the procedure is automated by the test software.

Approximate Test Time

1 to 2 hours depending upon the number of test temperatures needed.

Basic Procedure

• Heat the asphalt binder from which the test specimens are to be selected until the binder is sufficiently fluid to pour the test specimens (Video 1).

Cold asphalt binder can develop reversible molecular associations that cause it to stiffen (called “steric hardening”). Without heating, steric hardening can result in overestimating the complex modulus by as much as 50 percent (AASHTO, 2000c [1] ).

• Select the testing temperature according to the asphalt binder grade or testing schedule. Heat the DSR to the test temperature. This preheats the the upper and lower plates (Figure 9 and Figure 10), which allows the specimen to adhere to them.

• Place the asphalt binder sample between the test plates (Figure 11).

• Move the test plates together until the gap between them equals the test gap plus 0.002 inches (0.05 mm).
• Trim the specimen around the edge of the test plates using a heated trimming tool.

The calculated complex modulus (G*) is proportional to the fourth power of the asphalt binder specimen radius, therefore careful trimming will insure more reliable measurements (AASHTO, 2000c [1] ).

• Move the test plates together to the desired testing gap. This creates a slight bulge in the asphalt binder specimen’s perimeter.
• Bring the specimen to the test temperature. Start the test only after the specimen has been at the desired temperature for at least 10 minutes.
• The DSR software determines a target torque at which to rotate the upper plate based on the material being tested (e.g., unaged binder, RTFO residue or PAV residue). This torque is chosen to ensure that measurements are within the specimen’s region of linear behavior.
• The DSR conditions the specimen for 10 cycles at a frequency of 10 rad/sec (1.59 Hz).

Movement of the test plates at 10rad/sec is so small that you should not be able to easily see it. If movement is obvious, the bond between the asphalt binder sample and the test plates may have broken.

• The DSR takes test measurements over the next 10 cycles and then the software reduces the data to produce a value for complex modulus (G*) and phase angle (δ) (Figure 12).

Testing should be done as quickly as possible to minimize the effect of steric hardening that occurs during the test. Steric hardening can cause an increase in complex modulus (G*) if the specimen is kept in the DSR for a prolonged period of time.

When testing at multiple temperatures, all testing should be completed within four hours (AASHTO, 2000c [1] ).

Parameters Measured

• Complex modulus (G*)
• Phase angle (δ)

Specifications

Table 1: Performance Graded Asphalt Binder DSR specifications

Typical Values

The complex modulus (G*) can range from about 0.07 to 0.87 psi (500 to 6000 Pa), while the phase angle (δ) can ranges from about 50 to 90°. A δ of 90° is essentially complete viscous behavior. Polymer-modified asphalt binders generally exhibit a higher G* and a lower δ. This means they are, in general, a bit stiffer and more elastic than unmodified asphalt cements.

Calculations (see Interactive equation)

DSR software performs the necessary calculations automatically. The DSR software uses the following equations:

τ max = maximum applied stress

γ max = maximum resultant strain

T = maximum applied torque

r = specimen radius (either 4 or 12.5 mm)

θ = deflection (rotation) angle (in radians)

h = specimen height (either 1 or 2 mm)

δ = time lag between occurrence of τ max and γ max (see Figure 5)

The phase angle cannot be less than 0° or greater than 90°. The time lag can be measured in seconds and then converted to an angular measurement by dividing it by the oscillation frequency and then multiplying by 360° (or 2π radians).

• American Association of State Highway and Transportation Officials (AASHTO). (2000c).  AASHTO Provisional Standards, April 2000 Edition . American Association of State Highway and Transportation Officials. Washington, D.C. ↵

I. INTRODUCTION

Ii. experiments, a. experiment summary, b. cross-platform comparisons, iii. elastic constants, iv. results and discussion, a. shear modulus pressure scaling of strength, b. alternative mechanisms, v. conclusions, acknowledgments, data availability, experimental evaluation of shear modulus scaling of dynamic strength at extreme pressures.

• Split-Screen
• Article contents
• Figures & tables
• Supplementary Data
• Peer Review
• Open the PDF for in another window
• Reprints and Permissions
• Cite Icon Cite
• Search Site

J. L. Brown , M. B. Prime , N. R. Barton , D. J. Luscher , L. Burakovsky , D. Orlikowski; Experimental evaluation of shear modulus scaling of dynamic strength at extreme pressures. J. Appl. Phys. 28 July 2020; 128 (4): 045901. https://doi.org/10.1063/5.0012069

Download citation file:

• Ris (Zotero)
• Reference Manager

Recent progress in the development of dynamic strength experimental platforms is allowing for unprecedented insight into the assumptions used to construct constitutive models operating in extreme conditions. In this work, we make a quantitative assessment of how tantalum strength scales with its shear modulus to pressures of hundreds of gigapascals through a cross-platform examination of three dynamic strength experiments. Specifically, we make use of Split–Hopkinson pressure bar and Richtmyer–Meshkov instability experiments to assess the low-pressure strain and strain rate dependence. Concurrent examination of magnetically driven ramp-release experiments up to pressures of 350 GPa allows us to examine the pressure dependence. Using a modern description of the shear modulus, validated against both ab initio theory and experimental measurements, we then assess how the experimentally measured pressure dependence scales with shear modulus. We find that the common assumption of scaling strength linearly with the shear modulus is too soft at high pressures and offer discussion as to how descriptions of slip mediated plasticity could result in an alternative scaling that is consistent with the data.

Understanding the constitutive response of metals compressed to high energy density (HED) conditions, typically defined as pressures > 100 GPa ⁠ , is required for modeling a wide range of applications including high velocity impacts, planetary formation and interiors, and inertial confinement fusion implosions. 1 A key component to predictive models under HED conditions is a description of how the material strength scales with pressure. Early modeling efforts such as those of Steinberg et al. 2 utilized first order expansions to describe the pressure scaling of the shear modulus, G ⁠ , and strength, Y ⁠ . While it was recognized the derivatives of G and Y with respect to pressure could be different, most low pressure experiments examining temperature dependence suggest the local forces or stresses that must be overcome to drive dislocations past obstacles are linearly proportional to the shear modulus 3 such that

Elastic constants such as the shear modulus are significantly easier to determine through theory than the strength, so Eq. (1) has historically been a convenient way to describe the pressure scaling of the strength. Given the absence of any data to the contrary, similar approximations have been made in more modern strength models. The Preston–Tonks–Wallace (PTW) model, 4 for example, describes the strength as a function of strain, ϵ ⁠ , strain rate, ϵ ˙ ⁠ , pressure, P ⁠ , and temperature, T ⁠ , with an overall scaling by the shear modulus: Y = G ( P , T ) f ( ϵ , ϵ ˙ , T ) ⁠ . However, it is worth noting that f has a weak additional dependence on pressure through the calculation of a dimensionless strain rate variable. Similarly, in dislocation-based plasticity models such as that of Barton et al. , 5 the hardening contribution to the strength is generally described by the Taylor equation, Y ∝ G ϱ ⁠ , where ϱ is the dislocation density.

Recent progress in developing strength platforms on HED facilities such as the Z-machine 6 and the National Ignition Facility (NIF) 1 is allowing unprecedented access into the high pressure strength of materials. In this article, we explore the linear scaling approximation in Eq. (1) for the body-centered cubic (BCC) metal tantalum by examining experimental data from a range of platforms and comparing results to theoretical descriptions of the elastic constants. In Sec. II , we summarize previously published experimental results across multiple platforms and then use these data to estimate the strain, strain rate, temperature, and pressure dependencies. Section III examines first-principles calculations and an analytic model to estimate the high pressure behavior of the shear modulus. Sound speed measurements under both shockless (low temperature) and shock (high temperature) loading conditions validate the shear modulus model. Section IV describes how the high pressure strength data scale with the shear modulus and offer discussion on several physical mechanisms, which could influence the interpretation. Our conclusions are given in Sec. V .

In a coordinated effort between Los Alamos, Livermore, and Sandia national laboratories, a series of experiments were conducted to understand the response of tantalum (Ta) to a wide range of loading conditions. To facilitate cross-platform comparisons and eliminate uncertainties about the influence of microstructure, all samples were machined from a common lot of material. The initial density was measured to be ρ 0 = 16.65 ± 0.08 g / c m 3 ⁠ , while the ambient elastic longitudinal and shear wavespeeds were c L = 4.19 ± 0.03 and c S = 2.07 ± 0.02 km / s ⁠ , respectively. Based on these acoustic measurements, the shear modulus at room temperature and pressure was found to be G 0 = 71.3 ± 1.2 GPa ⁠ . Further microstructural characterization and quasi-static (QS) mechanical testing was reported previously by Buchheit et al. 7 In this section, we summarize previously reported data from three dynamic experimental platforms using this common lot of Ta and then make comparisons across the conditions accessed.

Table I summarizes a set of experiments spanning four decades in loading rates. First are Split–Hopkinson Pressure Bar (SHPB) experiments, which provide a stress–strain response at moderate ( ⁠ 10 3 s − 1 ⁠ ) strain rates under adiabatic compression near ambient pressures. The experiments listed in Table I are analogous to those described previously by Chen and Gray, 8 but were conducted using the common lot of Ta. Data from the highest rate room temperature experiment are shown in the inset in Fig. 1 .

Strain rate sensitivity of Ta at ambient pressure and temperature. Bayesian regression (linear in log–log space) between the SHPB and RMI experiments was used to estimate uncertainties at the characteristic Z strain rate. The inset shows how the SHPB data point was determined: the mean is calculated at a strain of 0.2 and the uncertainty (shaded region) accommodates strain hardening to strains of 0.4 ⁠ .

Summary of the nominal deformation conditions for each experiment.

Moving to significantly higher rates ( ⁠ 10 7 s − 1 ⁠ ), quantitative strength estimates have recently been extracted using the Richtmyer–Meshkov Instability (RMI) platform. 9 In these experiments, a gas gun produces a symmetric plate impact loading condition. Sinusoidal perturbations are machined in the rear free surface of the target plate, resulting in a Richtmyer–Meshkov instability upon shock breakout. The growth of the instability is measured using velocimetry and numerical simulations enable the estimation of material strength. Shots with shock magnitudes of 20 ⁠ , 30 ⁠ , and 34 GPa were conducted resulting in extreme strain rates and minimal heating. While the initial shock stresses are finite, the instability growth occurs at the free surface, so the strength near ambient pressure is believed to dominate the experimental response of interest.

The final set of experiments considered here is the ramp-release Z-machine results, which will be used to examine the high pressure strength scaling. Data were originally reported to 250 GPa ⁠ , 6 but a reanalysis along with new data extending to 350 GPa 9,10 will be used here. In these experiments, magnetic compression is used to shocklessly (ramp) compress the Ta at rates between those in the SHPB and RMI experiments. These experiments provide a continuous measurement of the bulk wavespeed (i.e., bulk modulus) up to a specified peak pressure and then the loading reversal is used to estimate the shear modulus and strength near the peak state. 11 Thus, the loading rates ( ⁠ ∼ 5 × 10 5 s − 1 ⁠ ) are comparable across different Z experiments, but the peak pressures, temperatures, and strains are different.

In order to evaluate the pressure scaling of the strength observed in the Z experiments, we must first estimate a Y 0 = Y ( ϵ ˙ ∼ 10 5 s − 1 , P = 0 ) ⁠ , which accounts for the proper strain and strain rate hardening. To do this, we utilize the SHPB and RMI data. Since these data are taken at approximately zero pressure and they bound the Z strain rate, interpolation can be performed to estimate Y 0 at the Z strain rates. The remaining complication is in accounting for the strain hardening. Fortunately, as detailed in Table I , there is significant overlap in the plastic strains across all three platforms.

The SHPB provides a continuous measure of stress as a function of strain, while for the RMI and Z shots, there is a characteristic strain associated with each measurement. In the Z experiments, the average strain varies from 0.2 at a peak pressure of 60 GPa to 0.45 at a peak of 350 GPa ⁠ . By examining the measured SHPB curve, we can assign a representative value and uncertainty, which incorporates the observed strain hardening over the range of strain accumulated in the Z experiments. The resulting estimate is 0.68 ± 0.06 GPa and, as shown by the shaded region in the inset in Fig. 1 , encompasses the strain hardening over the entirety of that measurement. While we do not make use of the data, the result from a quasi-static (QS) compression test is also shown in Fig. 1 . The assigned uncertainty of 0.06 GPa also captures the hardening observed over strains of 0.2–0.4 at QS loading rates. So, up to ∼ 10 3 s − 1 ⁠ , we are capturing the rate dependence of the strain hardening through the uncertainty estimate. In the case of the reported RMI values, the strength is an average extracted over perturbations of different sizes and, subsequently, incorporates strains ranging from 0.2 to 1.2. Thus, the error bars shown in Fig. 1 naturally incorporate strain hardening over the range of strains sampled in the Z experiments. While the RMI experiments span a larger range of strains that could bias the strength to higher values, we do not believe this is a significant effect; further discussion is given by Prime et al. 9 and in Sec. IV B .

Assuming interpolation between SHPB and RMI properly captures the strain and rate hardening, the effective Y 0 for the Z experiments can be estimated through the use of a simple scaling law. If the dominant strength contribution is from the Taylor hardening term and the dislocation density is near saturation, the result is power-law relationship between strength and rate: 4,12   log 10 ⁡ ( Y ) = a + b log 10 ⁡ ( ϵ ˙ ) ⁠ , where a and b are fitting parameters. Bayesian regression of the SHPB and RMI data over this model gives estimates of the probability distributions for a and b ⁠ , which is represented graphically by the shaded region in Fig. 1 . Subsequent evaluation of the regression model at the characteristic Z strain rates gives the estimate Y 0 = 1.02 ± 0.06 GPa ⁠ . As discussed by Prime et al. , 9 there are indications that the high rates associated with the RMI data are close to the transition to a dislocation drag dominated regime. This mechanism is discussed in more detail in Sec. IV B , but one consequence is the log-linear approximation may not hold across the regression range and the estimated Y 0 could be biased high. However, as will be shown, this potential biasing is conservative with respect to our conclusion about the shear modulus scaling, so the log-linear approximation does not alter our overall interpretation.

The thermodynamic variable we have not considered to this point is temperature, estimates of which are given in Table I . In the SHPB, a temperature increase comes from the amount of plastic work converted to heat. 8 Assuming a conversion of 100 % ⁠ , the temperature rise is only ∼ 100 K ⁠ . In addition to the plastic work heat, the Z and RMI experiments have additional thermal contributions from the loading path. While the RMI experiments include additional heating associated with shock dissipation, the Z experiments reach higher temperatures due to the larger amount of compression.

A model for estimating the thermal softening of the shear modulus is given by Eq. (3) in Sec. III . Jumping ahead, we can use this model to estimate the corresponding change in the strength assuming Y / G is constant; as discussed previously, this is a reasonable leading order approximation for the ambient pressures associated with Fig. 1 . The thermal model utilizes the temperature normalized by the melt temperature, values of which are also given in Table I . Propagating these ratios through Eq. (3) , we estimate thermal softening of the SHPB, RMI, and Z experiments to be ∼ 2 % − 3 % ⁠ , 3 % − 7 % ⁠ , and 2 % ⁠ , respectively. As these estimates are within the uncertainties shown in Fig. 1 , we do not believe temperature plays a significant role in the interpretation of strength in these experiments.

Evaluation of how strength scales with shear modulus requires a model for the shear modulus. Here, we examine a recent parameterization of an analytic model and compare to first-principles calculations to establish conservative variations on the model. The uncertainties are then validated against longitudinal sound speeds measured in shock experiments, assuming that the equation of state (EOS) is relatively well-known.

We use the model of Burakovsky–Greeff–Preston (BGP). 13 This model is constructed to capture behavior over the range of pressures and temperatures under consideration here, with the cold shear modulus given analytically by

The method for determining the parameters for the model, along with values for Ta can be found in Burakovsky et al. ; 14 reported values are G ( ρ 0 , 0 ) = 72.2 GPa ⁠ , ρ 0 = 16.74 g / c m 3 ⁠ , γ 1 = 3.0567 ⁠ , γ 2 = 1.0 ⁠ , and q = − 6.6 ⁠ .

The temperature dependence of the BGP model is given by

where β ( ρ ) = 0.23 + 0.21 ( ρ − ρ 0 ) 1.09 ⁠ . For the sake of brevity and because we are primarily concerned with low temperatures, we are not summarizing T m ( ρ ) ⁠ . However, the T m model form bears similarities with Eq. (2) and parameters can be found in Ref. 14 .

While the BGP model was partially fit using ab initio methods, it is useful to examine an independent set of ab initio calculations to assess variability in the theoretical approaches. Specifically, we use the first-principles-based work of Orlikowski et al. , 15 in which density functional theory (DFT) calculations combined with the quasiharmonic approximation for temperature effects were used to assess moduli over the full range of pressures and temperatures of interest here.

A comparison between the BGP model and first-principles calculations of the Ta shear modulus is given in Fig. 2 . The difference between the BGP and first-principles 300 K isotherms is ∼ 20 % ⁠ . As such, we assign a 20 % standard deviation (represented by the shaded region in Fig. 2 ) to the BGP model as a way to study variability in the BGP model. We emphasize that this is not a rigorous estimate of the uncertainty in the BGP model, rather it is a convenient way to propagate variability in different ab initio methods through to our pressure scaling conclusions.

Pressure dependence of the Ta shear modulus for two fundamental curves: the 300 K isotherm and the Hugoniot. Uncertainties of 20 % in the analytic BGP model are represented by the shaded regions and are assigned to incorporate the first-principles-based calculations.

Traditionally, it has proven extremely difficult to make high-precision measurements of the shear modulus at high pressures. The Z experiments provide measurements of the bulk and longitudinal wavespeeds at peak compression, but the subsequent shear modulus estimates have large uncertainties. As shown in Fig. 2 , the Z measurements are consistent with the BGP model (assuming the difference between the isotherm and quasi-isentrope are relatively small), but are not precise enough to validate the theory. An alternative type of experiment, which can provide a stronger constraint, are longitudinal sound speed measurements in shock-compressed material. Ta, in particular, has received a great deal of interest with experiments having been performed by multiple researchers. 16–20 However, to connect to these experiments, we require two additional pieces of information: the bulk modulus, B ⁠ , and the effect of temperature on the elastic constants ( ⁠ B and G ⁠ ) due to the shock heating.

The thermal description for the BGP model is given in Eq. (3) , while the first-principles-based method directly provides values along the Hugoniot. 15 The resulting shear moduli along the Hugoniot are given in Fig. 2 . As shown, there is significant softening of the shear modulus as Ta approaches shock melting ( ⁠ ∼ 300 GPa ⁠ ), and the BGP model and first-principles descriptions of this softening are in good agreement.

The longitudinal sound speed, c L ⁠ , is a combination of the bulk and shear response,

Thus, we require a model for the bulk modulus along the Hugoniot to tie to the sound speed measurements in shocked material. Toward this end, we use the Sesame 93524 model, 21 which has been validated against a variety of experimental and theoretical data. In Fig. 3 , we provide additional validation of the Sesame 93524 description for B against the Z experiments and first-principles-based calculations. The Z experiments provide a continuous measurement of the speed of sound/bulk modulus along a quasi-isentropic loading path, 6 which is in excellent agreement with both the first-principles and Sesame 93524 calculations of the isentrope to within the Z errors of ∼ 3 % ⁠ . The spike in Z measurement near zero pressure is a consequence of the initially elastic deformation, so this peak is representative of the longitudinal modulus before the rest of the measurement transitions to reflect the bulk response. Examining the Hugoniots, there is also good agreement between the first-principles calculations and the Sesame 93524 model. The thermal effects, diagnosed as the difference between the isentrope and Hugoniot, are shown to be significantly less on B than G ⁠ . As such, we have high confidence in the Sesame 93524 description for B ⁠ , making the sound speed measurements in shocked material a strong constraint on G ⁠ .

Pressure dependence of the Ta bulk modulus. The Sesame 93524 EOS is validated at low temperatures with the Z experiments (shaded region is the standard error) and at high temperatures with first-principles-based calculations. Given the temperature sensitivity and quality of data available (only a subset is shown here), there is significantly higher confidence in the bulk than shear modulus.

The BGP model for G is combined with the Sesame 93524 description of B through Eq. (4) to produce the model curve shown in Fig. 4 . Experimental data, 16–20 utilizing a variety of measurement techniques, are also shown. The BGP model with the associated 20 % uncertainties gives a good representation of the data. The notable exception is the data of Akin et al. 20 at pressures < 100 GPa ⁠ , but the source of this discrepancy is unclear. We further note the 1 σ uncertainties of 20 % capture nearly all of the data (treating the low pressure 20 data as an outlier), which suggests that the model uncertainties are overly conservative. Again, the BGP uncertainty is meant as an illustrative way to represent potential variability, so a conservative result within this validation is both desirable and expected.

Longitudinal speeds along the Ta Hugoniot have been measured in a variety of independent studies. 16–20 The models and first-principles-based points shown in Figs. 2 and 3 are combined through Eq. (4) to provide validation of the BGP representation of the shear modulus.

Sections II and III have provided the background to evaluate the linear pressure scaling assumption in Eq. (1) as it applies to tantalum. Here, we show that the linear scaling is not consistent with the data and offer an alternative scaling which is consistent. We then assess a variety of alternative physical mechanisms, which could result in an apparent non-linear scaling.

The BGP model for the shear modulus is described and validated in Sec. III and is shown in Fig. 2 . Propagating the isothermal shear modulus through Eq. (1) , along with the estimated uncertainty in Y 0 determined in Fig. 1 , gives the linear scaling result shown in Fig. 5 . As illustrated, this scaling is consistent with the Z data up to ∼ 100 GPa ⁠ , but diverges at higher pressures. Thus, even with a conservative estimate for uncertainty on the BGP model, the linear scaling by shear modulus approximation does not hold.

Scaling of the high pressure strength of tantalum with shear modulus. The estimated errors in Y 0 and G shown in Figs. 1 and 2 are propagated through the respective equations in the legend. Even with the conservative error estimates, the linear theory cannot explain the high pressure Z experiments. Instead, a non-linear theory is found to be consistent with the data. The dashed curve represents a hypothetical G ⁠ , which is fit such that the linear scaling approximation matches the data; as shown in Fig. 6 , this is not a plausible form for G ⁠ .

To rule out the possibility that an uncertainty in our shear modulus model could be leading to an erroneous interpretation, we fit a shear modulus, G hypothetical ⁠ , such that Eq. (1) reproduces the Z data. Thus, the dashed line in Fig. 5 is constructed to fit the data. The corresponding hypothetical shear modulus is compared with the validation data in Fig. 6 . Since G hypothetical is fit to the Z strength data it is interpreted to be representative of the Z quasi-isentropic loading path. As such, direct comparisons with the Z estimates of the shear modulus from Fig. 2 can be made, revealing poor agreement. More convincing is the translation to the sound speed data from the shock compression experiments summarized in Fig. 4 . To make this comparison, G hypothetical is translated to the Hugoniot using the same thermal models (a combination of the BGP temperature dependence in Eq. (3) and the Sesame 93524 Hugoniot) used to generate the analytic model in Fig. 4 . As shown, there is an irreconcilable difference with the experimental measurements at higher pressures. By extension, this also makes G hypothetical inconsistent with the first-principles calculations. As such, we conclude that G hypothetical is not a plausible explanation to describe the discrepancy observed in the pressure scaling result in Fig. 5 .

Evaluation of the hypothetical shear modulus against experimental measurements. G hypothetical is constructed such that a linear scaling of the strength provides a fit to the Z strength data (shown in Fig. 5 ). This G hypothetical is found to be significantly higher than shear modulus estimates from the Z experiments. Similarly, translating to a sound speed in the shocked state, the hypothetical solution is unreasonably high when compared to the experimental data (shown in Fig. 4 ).

In making these pressure scaling arguments, it is assumed that the Z experiments have been interpreted correctly to produce a meaningful inference of strength, so an independent dataset would be valuable in validating these results. Providentially, significant effort has been put into performing a completely different type of dynamic high pressure strength experiment at the National Ignition Facility (NIF), which provides an opportunity for cross-platform validation of the pressure scaling. These experiments monitor the growth of a Rayleigh–Taylor (RT) instability and have been performed to comparable pressures as the Z experiments (350 GPa), albeit at larger plastic strains and strain rates. Initial strength estimates of ∼ 10 GPa at a pressure of ∼ 350  GPa and a strain rate of ∼ 10 7 s − 1 were reported for the RT experiments. 1 Given the rate dependence described in Sec. II B , Y 0 for the RT experiments can be approximated by the RMI data, which have a mean of 1.3 GPa. So, to first order, the RT experiments would be expected to have a strength ∼ 30 % higher than the Z experiments at the same pressure. Given the uncertainties on the Z data are about this magnitude, the comparable values of strength between the two platforms agree to within the error bars. Thus, initial indications suggest the NIF RT experiments are consistent with the non-linear pressure scaling observed in the Z data. However, extracting a characteristic value from the RT experiments is extremely challenging, and, as such, we are in the process of a more thorough evaluation of the RT data to better assess this claim.

There are a variety of reasons the pressure scaling of the strength could be larger than the simple linear shear modulus scaling. Plasticity is typically accommodated by the motion of dislocations through the lattice; thus, any effects on the formation, structure, and motion of dislocations will have an effect on the strength. In general this behavior is significantly path dependent, because the behavior of individual dislocations depends upon overall dislocation structure, and dislocation structure in turn evolves throughout the deformation history. In the following paragraphs, we provide speculative examples of phenomena that may not be captured by the simple shear modulus scaling. However, we are not asserting that any one of these is the actual explanation of the discrepancy.

As a first example, we examine a traditional view of plastic slip, which comprises a summation of the Peierls stress and an additional strength term due to dislocation interactions. Throughout this article, we refer to this latter term as Taylor hardening. Because higher strain rates cause a higher saturation dislocation density, it is plausible the Taylor hardening becomes more significant at higher strain rates. Thus, it is possible that the Taylor hardening term has pressure effects that go beyond the traditional shear modulus scaling; for example, pressure could influence the dislocation kinetics in a manner that produces increased saturation dislocation density.

It is also possible that pressure would have an effect on the dislocation core structure and thus mobility even of pure screw dislocations. Correspondingly, the dislocation response could depend on the local stress field in a more complex manner, resulting, in some cases, in an additional pressure dependence of the strength. Furthermore, because the Peierls barrier for screw dislocation motion in tantalum is large, the motion of screw dislocations is accommodated by the formation of pairs of “kinks” in the otherwise straight screw dislocation. These kinks are of mixed-character (that is, include some edge dislocation) and are of higher mobility than the pure screw dislocation. Thus, any mechanism with an effect on the motion of edge dislocations may also affect the motion of screw dislocations.

Increasing pressure may suppress kink-pair nucleation due to the dilatational nature of the mixed-character kink pairs. Couch and Swartz, 22 for example, computed the dilatancy associated with dislocation kink pairs in BCC lattices and related this to an activation volume associated with vacancy-type point defects. Because of the dilatation of the lattice strain field associated with dislocations, there is a relationship between the defect energy associated with dislocations, the applied pressure, and the mobility of dislocations. 23–25 The generation of dislocations and their subsequent motion causes a change in the local crystal volume and, consequently, an increased amount of work must be applied to move the defect in the presence of a resisting pressure field. 23 Jung 23 performed an analysis to demonstrate that this mobility effect may result in the strength behaving as

For a linear variation of shear modulus with pressure, this relationship would exhibit an increase in strength with respect to pressure that is twice that exhibited by Eq. (1) . This factor of two represents the interaction of dilatational dislocation strain fields with pressure, producing an increased magnitude of the Peierls barrier. Figure 5 shows the result of using the Eq. (5) scaling relationship with our estimated Y 0 and the BGP model, with the result being in excellent agreement with the higher pressure Z data. Thus, we offer Eq. (5) as a plausible expression for the pressure scaling of Ta to 350 GPa ⁠ . More generally, comparing the first order expansion of Eq. (1) with Eq. (5) suggests that the pressure dependence of the strength may scale as Y ≈ Y 0 [ 1 + n P G 0 ∂ G ∂ P ] ⁠ , where n is a material-dependent parameter and our results indicate n ≈ 2 for Ta.

In Sec. IV A , we gave examples of why slip-based plasticity could have a non-linear scaling with pressure. Here, we examine other physical mechanisms that could explain the apparent discrepancy of the simple linear scaling shown in Fig. 5 and present arguments as to why we believe these mechanisms are a less likely explanation. The fact that the experiments conducted on Z are shockless (quasi-isentropic ramp compression) and induce strain rates on the order of 5 × 10 5 s − 1 ⁠ , immediately makes many of the proposed mechanisms highly unlikely. However, in the interest of broader discussion, it is valuable to consider how these other effects associated with the path-dependence of the material response can be masked when considering only the pressure and strain rate.

1. Strain hardening

First, we revisit the strain hardening discussion in Sec. II B in the context of pressure scaling. In order to interpret the observed dependence of strength on pressure in the Z experiments, it is necessary to understand other dependencies that are not depicted by the pressure axis alone. Consider that the Z experiments probe a 1D uniaxial deformation such that the increases in density (and pressure) are nearly directly associated with corresponding increases in the plastic strain. As higher pressures are probed, the plastic strain experienced by the material is also increased. Thus, there is a convolution of pressure and strain hardening effects in these data. Because of this correlation between pressure and plastic strain, discrepancy between the linear model and the data ( Fig. 5 ) could be associated with either the strain hardening behavior, the pressure dependence of the strength, or both. As such, the method described in Sec. II B to determine uncertainties on the representative Y 0 incorporated the strain hardening behavior over the range of plastic strains accessed in the Z experiments. As described in Table I , the peak plastic strain varies from 0.2 to 0.45 in the lowest to highest pressure Z experiments, respectively. The mean Y 0 for the SHPB bar data was taken at a strain of 0.2 ⁠ , but uncertainties reflect the strain hardening spanning the range of the data, up to strains of 0.4 ⁠ . Similarly, the RMI data points reflect a range of strains encompassing the Z data. Therefore, the interpolation to Y 0 for the Z strain rate encompasses the strain hardening of the two bounding experiments. We further note that recent advances in miniature Hopkinson bars have provided Ta measurements at strain rates comparable to those in the Z experiments. 26 In these experiments, the strain hardening is slight and comparable to the trend past strains of 0.1 shown in Fig. 1 . Thus, a variety of experiments suggest a small degree of strain hardening over the strains sampled in the Z data and the amount of hardening has been quantified and incorporated within the estimate used for Y 0 ⁠ . As such, if traditional assumptions about strain hardening hold, then we contend that strain hardening has been accounted for in the interpretation of Fig. 5 . As noted in Sec. IV A , however, it is possible that strain hardening has unknown pressure dependencies and further work would be needed to evaluate this possibility.

2. Non-Schmid behavior

The driving force for plastic slip in FCC metals is classically understood to be the shear stress projected onto each particular crystallographic slip plane, i.e., Schmid stress. However, in BCC metals including Ta, it is often reported that there is an additional resistance to (or enhancement of) slip on particular crystallographic planes associated with the non-Schmid components of the stress. 27–30 These non-Schmid effects manifest as a dependence of the apparent slip resistance on stress components that do not contribute to the Peach–Koehler force acting on a particular dislocation. 31 There are various phenomenological explanations of this apparent behavior; for example, the dislocation core with Burgers vector [111] spreads onto several intersecting slip planes that also contain the [111] vector, and non-Schmid stress components modify this dislocation core structure to affect the energy barrier for dislocation glide along the slip plane. 32,33 Alternately, the non-Schmid stress components may affect the nucleation of kink pairs; thus, as discussed in Sec. IV A , modifying the apparent mobility of screw dislocations. 34 There are a variety of theories to incorporate this behavior into continuum models which add additional non-Schmid stress components to the driving force for dislocation-accommodated crystallographic slip. 35–39 Most of these models are constrained such that hydrostatic stress does not influence the yield behavior. However, Yalcinkaya et al. 40 and Cho et al. 39 used a five-term non-Schmid stress yield criterion, which does include a pressure dependence. While non-Schmid behavior can manifest in asymmetric strength behavior under tension vs compression, the effect is more closely associated with dislocation core mobility under shearing stresses and typically not with the applied pressure. In this paper, the focus is on the relationship between observed strength and pressure. Furthermore, the non-Schmid behavior and associated parameters are known to depend upon temperature 30,37,41–43 with stronger effects at cold temperatures (4–300 K) and only minimal evidence of the behavior at the temperatures relevant to the high-pressure experiments evaluated here.

3. Phonon drag

At sufficiently high strain rates (well above the Z rates), dislocation velocities are governed by the interaction of dislocations with lattice phonons. This interaction is of a viscous nature, which is characterized by the phonon drag viscosity. Many models of the phonon drag viscosity increase with temperature and, nominally, also vary with density via a term that is largely related to the longitudinal wave speed. The details of this relationship could be strongly influenced by the relationship between second- and third-order elastic constants at varying densities. 44,45 For the purposes of discussion, it is plausible that the strength in the viscous drag regime does not scale linearly with shear modulus, although such effects have received less attention in the literature.

As discussed previously, there is evidence the RMI data may be close to this transition to a phonon drag regime, so there is a possibility the fitting in Sec. II B biased the Y 0 in Fig. 5 to too high of a value. Fortunately, this does not change our interpretation of the shear modulus scaling relationship as a lower Y 0 would only push the linear scaling further from the data.

4. Homogeneous dislocation nucleation

Homogeneous nucleation of dislocations may be expected at shock pressures in excess of 50 GPa ⁠ . Since the Z data are shockless; however, this mechanism is also unlikely. Dislocation nucleation has been observed in some molecular dynamics (MD) simulations, 46–48 but is noted to be dependent on the specific potential used. Furthermore, recent MD calculations with sufficiently large volumes of material have shown that the material may tend to undergo dislocation multiplication or twinning rather than dislocation nucleation. 49 Overall, it is unlikely that dislocation nucleation plays a significant role in these ramp compression-release experiments.

5. Twinning

As with homogeneous dislocation nucleation, deformation by mechanical twinning under shockless compression is unlikely. Negligible twinning has been observed in single crystal tantalum shocked to ∼ 25 GPa in both recovery 50 and in situ x-ray diffraction experiments. 51 However, in both of these works, significant twinning was observed at higher pressures ( ⁠ > 50 GPa ⁠ ). 50,51 Interestingly, the work of Wehrenberg et al. 51 suggests a transition back to dislocation-slip-dominated plasticity at pressures beyond 150 GPa ⁠ . In additional recovery work on polycrystalline material, similar trends are observed; however, even modest amounts of initial cold work, and the corresponding increase in dislocation density, suppress most twinning. 52 The progression of twin evolution appears to be a complex phenomenon, with modeling work suggesting that, due to localized stress fields produced in the vicinity of the twins, once twins have formed, they may propagate even at macroscopic conditions under which they would not have formed. 53 Overall, whether or not mechanical twinning is a significant deformation mechanism in tantalum is expected to depend at least on prior processing (initial dislocation structure) and the magnitude (and even character 54 ) of any initial shock. Given the nature of the common lot of Ta material (which was cold-rolled 7 ) and the shockless deformation path of the Z experiments, twinning is not expected have significant bearing on the results examined here.

6. Phase transformation

The final potential mechanism of note is a pressure induced phase transformation, which would obviously invalidate any arguments about scaling from properties of the ambient phase. The hexagonal omega phase has been observed in shock-recovered Ta driven to 45 GPa with explosively driven flyer plates 55 and to 70 GPa with the Omega laser facility. 56 Correspondingly, some ab initio DFT calculations suggest a polymorphic phase transformation near the melt line, and this phase boundary is estimated to intersect the Hugoniot at ∼ 100 GPa ⁠ . 57 As with previously discussed mechanisms, the shockless loading of the Z experiments likely negates any concerns about this potential transition. The room temperature isotherm is a significantly better approximation to the Z loading path than the Hugoniot, particularly at high pressures. In this case, DFT calculations of the cold curve suggest that Ta remains in its ambient phase to 1000 GPa 15 and x-ray diffraction measurements in static diamond anvil cell experiments validate this prediction up to 300 GPa ⁠ . 58  

The high pressure strength scaling of Ta was assessed through a critical examination of several dynamic experimental platforms. We examined two low-pressure platforms, Split–Hopkinson pressure bar and Richtmyer–Meshkov instability, to isolate hardening effects and strain rate sensitivity. The measured response was then mapped to a third platform, Z ramp-release, allowing for the estimation of the zero-pressure strength at the characteristic Z strain rate. Using the high-pressure Z measurements (up to 350 GPa), we were then able to assess various assumptions about how the strength scales with shear modulus. We find that the most common approximation of scaling the strength linearly with the shear modulus does not hold beyond pressures of ∼ 100  GPa. Through an examination of other potential mechanisms that could explain this apparent discrepancy, such as twinning, homogeneous dislocation nucleation, or phase transformation, we conclude that the most plausible deformation mechanism is slip mediated plasticity, which suggests that a non-linear dependence of the strength with shear modulus is required to model the experimental data.

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under Contract No. DE-NA-0003525. Parts of this work were supported by the U.S. Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 189233218CNA000001). Parts of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA2734. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

We would like to thank the rest of the tri-lab team for the discussions that inspired this work. We would also like to thank Dan Casem of the U.S. Army Research Laboratory, Aberdeen, for helpful discussions on the measured high-rate behavior of Ta.

The data that support the findings of this study are available within the article.

Citing articles via

Submit your article.

Sign up for alerts

Related Content

• Online ISSN 1089-7550
• Print ISSN 0021-8979
• For Researchers
• For Librarians
• For Advertisers
• Our Publishing Partners  
• Physics Today
• Conference Proceedings
• Special Topics

pubs.aip.org

• Privacy Policy
• Terms of Use

Connect with AIP Publishing

This feature is available to subscribers only.

Sign In or Create an Account

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

• View all journals
• Explore content
• About the journal
• Publish with us
• Sign up for alerts
• Published: 24 September 2024

Dynamic traction force measurements of migrating immune cells in 3D biopolymer matrices

• David Böhringer   ORCID: orcid.org/0000-0002-3880-4069 1 ,
• Mar Cóndor   ORCID: orcid.org/0000-0002-8656-7846 2 , 3 ,
• Lars Bischof 1 ,
• Tina Czerwinski 1 ,
• Niklas Gampl   ORCID: orcid.org/0009-0002-1392-0108 4 , 5 ,
• Phuong Anh Ngo 6 , 7 ,
• Andreas Bauer 1 ,
• Caroline Voskens 7 , 8 ,
• Rocío López-Posadas 6 , 7 ,
• Kristian Franze 4 , 5 , 9 ,
• Silvia Budday   ORCID: orcid.org/0000-0002-7072-8174 10 ,
• Christoph Mark   ORCID: orcid.org/0000-0002-8612-6469 1 ,
• Ben Fabry   ORCID: orcid.org/0000-0003-1737-0465 1 &
• Richard Gerum   ORCID: orcid.org/0000-0001-5893-2650 1 , 11  

Nature Physics ( 2024 ) Cite this article

1 Altmetric

Metrics details

• Biological physics
• Cellular motility
• Computational biophysics

Immune cells, such as natural killer cells, migrate with high speeds of several micrometres per minute through dense tissue. However, the magnitude of the traction forces during this migration is unknown. Here we present a method to measure dynamic traction forces of fast migrating cells in biopolymer matrices from the observed matrix deformations. Our method accounts for the mechanical nonlinearity of the three-dimensional tissue matrix and can be applied to time series of confocal or bright-field image stacks. It allows for precise force reconstruction over a wide range of force magnitudes and object sizes—even when the imaged volume captures only a small part of the matrix deformation field. We demonstrate the broad applicability of our method by measuring forces from around 1 nN for axon growth cones up to around 10 μN for mouse intestinal organoids. We find that natural killer cells show bursts of large traction forces around 50 nN that increase with matrix stiffness. These force bursts are driven by myosin II contractility, mediated by integrin β1 adhesions, focal adhesion kinase and Rho-kinase activity, and occur predominantly when the cells migrate through narrow matrix pores.

This is a preview of subscription content, access via your institution

Access options

Access Nature and 54 other Nature Portfolio journals

Get Nature+, our best-value online-access subscription

24,99 € / 30 days

cancel any time

Subscribe to this journal

Receive 12 print issues and online access

251,40 € per year

only 20,95 € per issue

Buy this article

• Purchase on SpringerLink
• Instant access to full article PDF

Prices may be subject to local taxes which are calculated during checkout

Similar content being viewed by others

Microparticle traction force microscopy reveals subcellular force exertion patterns in immune cell–target interactions

Slow integrin-dependent migration organizes networks of tissue-resident mast cells

Detecting long-range interactions between migrating cells

Data availability.

Source data from this study are are provided with this paper and as part of the Saenopy GitHub repository 50 . Image raw data are available upon request from the corresponding author.

Code availability

The software (Saenopy) and all dependent packages are available on GitHub as an open-source Python package with a graphical user interface 50 . Figures are created using the Python packages Pylustrator 66 and PyVista 67 using a colourblind-friendly colour palette from ref. 68 .

Dembo, M. & Wang, Y. L. Stresses at the cell-to-substrate interface during locomotion of fibroblasts. Biophys. J. 76 , 2307–2316 (1999).

Article   ADS   Google Scholar  

Han, S. J., Oak, Y., Groisman, A. & Danuser, G. Traction microscopy to identify force modulation in subresolution adhesions. Nat. Methods 12 , 653–656 (2015).

Article   Google Scholar  

Huang, Y. et al. Traction force microscopy with optimized regularization and automated Bayesian parameter selection for comparing cells. Sci. Rep. 9 , 539 (2019).

Butler, J. P., Tolic-Norrelykke, I. M., Fabry, B. & Fredberg, J. J. Traction fields, moments, and strain energy that cells exert on their surroundings. Am. J. Physiol. Cell Physiol. 282 , C595–C605 (2002).

Bergert, M. et al. Confocal reference free traction force microscopy. Nat. Commun. 7 , 12814 (2016).

Sabass, B., Gardel, M. L., Waterman, C. M. & Schwarz, U. S. High resolution traction force microscopy based on experimental and computational advances. Biophys. J. 94 , 207–220 (2008).

Bauer, A. et al. pyTFM: a tool for traction force and monolayer stress microscopy. PLoS Comput. Biol. 17 , e1008364 (2021).

Gómez-González, M., Latorre, E., Arroyo, M. & Trepat, X. Measuring mechanical stress in living tissues. Nat. Rev. Phys. 2 , 300–317 (2020).

Barrasa-Fano, J. et al. TFMLAB: a MATLAB toolbox for 4D traction force microscopy. SoftwareX 15 , 100723 (2021).

Legant, W. R. et al. Measurement of mechanical tractions exerted by cells in three-dimensional matrices. Nat. Methods 7 , 969–971 (2010).

Franck, C., Maskarinec, S. A., Tirrell, D. A. & Ravichandran, G. Three-dimensional traction force microscopy: a new tool for quantifying cell–matrix interactions. PLoS ONE 6 , e17833 (2011).

Steinwachs, J. et al. Three-dimensional force microscopy of cells in biopolymer networks. Nat. Methods 13 , 171–176 (2016).

Song, D. et al. Recovery of tractions exerted by single cells in three-dimensional nonlinear matrices. J. Biomech. Eng. 142 , 081012 (2020).

Dong, L. & Oberai, A. A. Recovery of cellular traction in three-dimensional nonlinear hyperelastic matrices. Comput. Methods Appl. Mech. Eng. 314 , 296–313 (2017).

Article   ADS   MathSciNet   Google Scholar  

Toyjanova, J. et al. High resolution, large deformation 3D traction force microscopy. PLoS ONE 9 , e90976 (2014).

Storm, C., Pastore, J. F., MacKintosh, F. C., Lubensky, T. C. & Janmey, P. A. Nonlinear elasticity in biological gels. Nature 435 , 188–191 (2005).

Munster, S. et al. Strain history dependence of the nonlinear stress response of fibrin and collagen networks. Proc. Natl Acad. Sci. USA 110 , 12197–12202 (2013).

Discher, D. E., Janmey, P. & Wang, Y. L. Tissue cells feel and respond to the stiffness of their substrate. Science 310 , 1139–1143 (2005).

Friedl, P., Noble, P. B. & Zänker, K. S. T lymphocyte locomotion in a three-dimensional collagen matrix. Expression and function of cell adhesion molecules. J. Immunol. 154 , 4973–4985 (1995).

Friedl, P., Entschladen, F., Conrad, C., Niggemann, B. & Zänker, K. S. CD4 + T lymphocytes migrating in three-dimensional collagen lattices lack focal adhesions and utilize β1 integrin-independent strategies for polarization, interaction with collagen fibers and locomotion. Eur. J. Immunol. 28 , 2331–2343 (1998).

3.0.CO;2-C" data-track-item_id="10.1002/(SICI)1521-4141(199808)28:08 3.0.CO;2-C" data-track-value="article reference" data-track-action="article reference" href="https://doi.org/10.1002%2F%28SICI%291521-4141%28199808%2928%3A08%3C2331%3A%3AAID-IMMU2331%3E3.0.CO%3B2-C" aria-label="Article reference 20" data-doi="10.1002/(SICI)1521-4141(199808)28:08 3.0.CO;2-C">Article   Google Scholar  

Yamada, K. M. & Sixt, M. Mechanisms of 3D cell migration. Nat. Rev. Mol. Cell Biol. 20 , 738–752 (2019).

Caillier, A., Oleksyn, D., Fowell, D. J., Miller, J. & Oakes, P. W. T cells use focal adhesions to pull themselves through confined environments. J. Cell Biol. 223 , e202310067 (2024).

Kaltenbach, L. & Lämmermann, T. Mast cell migration and organization in tissues depend on integrin–ECM interactions. Nat. Immunol. 24 , 899–900 (2023).

François, J. et al. The interplay between matrix deformation and the coordination of turning events governs directed neutrophil migration in 3D matrices. Sci. Adv. 7 , eabf3882 (2021).

Huse, M. Mechanical forces in the immune system. Nat. Rev. Immunol. 17 , 679–690 (2017).

Cóndor, M., Steinwachs, J., Mark, C., García-Aznar, J. M. & Fabry, B. Traction force microscopy in 3-dimensional extracellular matrix networks. Curr. Protoc. Cell Biol. 75 , 10.22.1–10.22.20 (2017).

Liberzon, A., Käufer, T., Bauer, A., Vennemann, P. & Zimmer, E. OpenPIV/openpiv-python: OpenPIV-Python v0.23.4. Zenodo https://doi.org/10.5281/zenodo.4409178 (2021).

Thievessen, I. et al. Vinculin is required for cell polarization, migration, and extracellular matrix remodeling in 3D collagen. FASEB J. 29 , 4555–4567 (2015).

Cóndor, M. et al. Breast cancer cells adapt contractile forces to overcome steric hindrance. Biophys. J. 116 , 1305–1312 (2019).

Wolf, K., Müller, R., Borgmann, S., Bröcker, E. B. & Friedl, P. Amoeboid shape change and contact guidance: T-lymphocyte crawling through fibrillar collagen is independent of matrix remodeling by MMPs and other proteases. Blood 102 , 3262–3269 (2003).

Czerwinski, T. et al. Immune cells employ intermittent integrin-mediated traction forces for 3D migration. Preprint at bioRxiv https://doi.org/10.1101/2023.04.20.537658 (2023).

Reinhart-King, C. A., Dembo, M. & Hammer, D. A. Cell–cell mechanical communication through compliant substrates. Biophys. J. 95 , 6044–6051 (2008).

Paszek, M. J. et al. Tensional homeostasis and the malignant phenotype. Cancer Cell 8 , 241–254 (2005).

Heidemann, S. R., Lamoureux, P. & Buxbaum, R. E. Growth cone behavior and production of traction force. J. Cell Biol. 111 , 1949–1957 (1990).

Betz, T., Koch, D., Lu, Y. B., Franze, K. & Käs, J. A. Growth cones as soft and weak force generators. Proc. Natl Acad. Sci. USA 108 , 13420–13425 (2011).

Koch, D., Rosoff, W. J., Jiang, J., Geller, H. M. & Urbach, J. S. Strength in the periphery: growth cone biomechanics and substrate rigidity response in peripheral and central nervous system neurons. Biophys. J. 102 , 452–460 (2012).

Polackwich, R. J., Koch, D., McAllister, R., Geller, H. M. & Urbach, J. S. Traction force and tension fluctuations in growing axons. Front. Cell. Neurosci. 9 , 417 (2015).

Santos, T. E. et al. Axon growth of CNS neurons in three dimensions is amoeboid and independent of adhesions. Cell Rep. 32 , 107907 (2020).

Singer, A. J. & Clark, R. A. Cutaneous wound healing. N. Engl. J. Med. 341 , 738–746 (1999).

Li, B. & Wang, J. H. Fibroblasts and myofibroblasts in wound healing: force generation and measurement. J. Tissue Viability 20 , 108–120 (2011).

Kang, N. Mechanotransduction in liver diseases. Semin. Liver Dis. 40 , 84–90 (2020).

Joyce, J. A. & Pollard, J. W. Microenvironmental regulation of metastasis. Nat. Rev. Cancer 9 , 239–252 (2009).

Calvo, F. et al. Mechanotransduction and YAP-dependent matrix remodelling is required for the generation and maintenance of cancer-associated fibroblasts. Nat. Cell Biol. 15 , 637–646 (2013).

Sato, T. et al. Long-term expansion of epithelial organoids from human colon, adenoma, adenocarcinoma, and Barrett’s epithelium. Gastroenterology 141 , 1762–1772 (2011).

Eiraku, M. et al. Self-organizing optic-cup morphogenesis in three-dimensional culture. Nature 472 , 51–58 (2011).

Lancaster, M. A. et al. Cerebral organoids model human brain development and microcephaly. Nature 501 , 373–379 (2013).

Gjorevski, N. et al. Tissue geometry drives deterministic organoid patterning. Science 375 , eaaw9021 (2022).

Pérez-González, C. et al. Mechanical compartmentalization of the intestinal organoid enables crypt folding and collective cell migration. Nat. Cell Biol. 23 , 745–757 (2021).

Poling, H. M. et al. Mechanically induced development and maturation of human intestinal organoids in vivo. Nat. Biomed. Eng. 2 , 429–442 (2018).

Gerum, R., Böhringer, D. & Cóndor, M. Saenopy Github Project. GitHub https://github.com/rgerum/saenopy (2020).

Bonferroni, C. E. Teoria Statistica Delle Classi e Calcolo Delle Probabilità (Florence R. Istituto Superiore di Scienze Economiche e Commerciali, 1936).

Steinwachs, J. Cellular Forces during Migration through Collagen Networks . PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg (2015).

Kah, D. et al. A low-cost uniaxial cell stretcher for six parallel wells. HardwareX 9 , e00162 (2021).

Faust, U. et al. Cyclic stress at mHz frequencies aligns fibroblasts in direction of zero strain. PLoS ONE 6 , e28963 (2011).

Koser, D. E. et al. Mechanosensing is critical for axon growth in the developing brain. Nat. Neurosci. 19 , 1592–1598 (2016).

Nieuwkoop, P. D. & Faber, J. Normal Table of Xenopus laevis (Daudin): A Systematical and Chronological Survey of the Development from the Fertilized Egg till the End of Metamorphosis (Garland, 1994).

Mark, C. et al. Collective forces of tumor spheroids in three-dimensional biopolymer networks. eLife 9 , e51912 (2020).

Martínez-Sánchez, L. D. C. et al. Epithelial RAC1-dependent cytoskeleton dynamics controls cell mechanics, cell shedding and barrier integrity in intestinal inflammation. Gut 72 , 275–294 (2023).

Otsu, N. A threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern. 9 , 62–66 (1979).

Tikhonov, A. N., Goncharsky, A. V., Stepanov, V. V. & Yagola, A. G. Numerical Methods for the Solution of Ill-Posed Problems (Springer Science & Business Media, 1995).

Tikhonov, A. N. On the stability of inverse problems. Proc. USSR Acad. Sci. 39 , 195–198 (1943).

MathSciNet   Google Scholar  

Huber, P. J. Robust Statistics Vol. 523 (John Wiley & Sons, 2004).

Kruskal, W. H. & Wallis, W. A. Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47 , 583–621 (1952).

Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17 , 261–272 (2020).

Harris, C. R. et al. Array programming with NumPy. Nature 585 , 357–362 (2020).

Gerum, R. Pylustrator: code generation for reproducible figures for publication. J. Open Source Softw. 5 , 1989 (2020).

Sullivan, C. & Kaszynski, A. PyVista: 3D plotting and mesh analysis through a streamlined interface for the Visualization Toolkit (VTK). J. Open Source Softw. 4 , 1450 (2019).

Wong, B. Points of view: color blindness. Nat. Methods 8 , 441 (2011).

Kroy, K. & Frey, E. Force-extension relation and plateau modulus for wormlike chains. Phys. Rev. Lett. 77 , 306–309 (1996).

Jansen, K. A. et al. The role of network architecture in collagen mechanics. Biophys. J. 114 , 2665–2678 (2018).

Lang, N. R. et al. Estimating the 3D pore size distribution of biopolymer networks from directionally biased data. Biophys. J. 105 , 1967–1975 (2013).

Krauss, P., Metzner, C., Lange, J., Lang, N. & Fabry, B. Parameter-free binarization and skeletonization of fiber networks from confocal image stacks. PLoS ONE 7 , e36575 (2012).

Sato, Y. et al. Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images. Med. Image Anal. 2 , 143–168 (1998).

Yen, J. C., Chang, F. J. & Chang, S. A new criterion for automatic multilevel thresholding. IEEE Trans. Image Process. 4 , 370–378 (1995).

Download references

Acknowledgements

This work was funded by the German Research Foundation (DFG; project 326998133 – TRR-SFB 225 – subprojects A01, B09 and C02; project 460333672 – CRC 1540 – subprojects A01, A05, B02 and C05; project 375876048 – TRR-SFB 241 – subprojects A07 and C04; project 461063481 – LO 2465/6-1; project 414058251 – SPP-1782, LO 2465/2-1), the National Institutes of Health (HL120839), the Emerging Fields Initiative of the University of Erlangen-Nuremberg and the Alexander von Humboldt Foundation (Humboldt Professorship, K.F.). We thank the ENB Biological Physics programme of the University Bayreuth for support, R. Henriques for the LaTeX layout, I. Thievessen for help with differential interference contrast imaging and R. Reimann for designing the Saenopy logo.

Author information

Authors and affiliations.

Department of Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany

David Böhringer, Lars Bischof, Tina Czerwinski, Andreas Bauer, Christoph Mark, Ben Fabry & Richard Gerum

Life Sciences Technology Department, Interuniversity Micro-Electronics Centre, Leuven, Belgium

Biomechanics Section, Department of Mechanical Engineering, KU Leuven, Leuven, Belgium

Institute of Medical Physics and Microtissue Engineering, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany

Niklas Gampl & Kristian Franze

Max-Planck-Zentrum für Physik und Medizin, Erlangen, Germany

Department of Medicine 1, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany

Phuong Anh Ngo & Rocío López-Posadas

Deutsches Zentrum Immuntherapie (DZI), Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany

Phuong Anh Ngo, Caroline Voskens & Rocío López-Posadas

Department of Dermatology, University Hospital Erlangen, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany

Caroline Voskens

Department of Physiology, Development and Neuroscience, University of Cambridge, Cambridge, UK

Kristian Franze

Institute of Continuum Mechanics and Biomechanics, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany

Silvia Budday

Department of Physics and Astronomy, York University, Toronto, Ontario, Canada

Richard Gerum

You can also search for this author in PubMed   Google Scholar

Contributions

Methodology: R.G., D.B., M.C., B.F. and C.M. Software: R.G., D.B., M.C. and A.B. Rheology: D.B., M.C., L.B., S.B. and N.G. Immune cell experiments: L.B., T.C., D.B. and C.V. Organoid experiments: P.A.N., D.B. and R.L.-P. Neuronal cell experiments: K.F. and N.G. Data analysis: D.B., L.B., C.M. and R.G. Writing: D.B., B.F., M.C., C.M. and R.G.

Corresponding author

Correspondence to Ben Fabry .

Ethics declarations

Competing interests.

The authors declare no competing interests.

Peer review

Peer review information.

Nature Physics thanks Hawa Racine Thiam, Karin Wang and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Extended data

Extended data fig. 1 nonlinear material model..

The nonlinear material model divides the mechanical response of individual fibers into a region where fiber stiffness ( ω ″ ) decreases exponentially with decreasing strain under compression (buckling), a region of constant fiber stiffness for small strains (straightening), and a region of exponentially increasing fiber stiffness for larger strains (stretching) 12 . Non-linear materials are modelled according to Eq. ( 1 ) using four material parameters: The linear stiffness k 0 , the buckling coefficient d 0 , the characteristic strain λ s , and the stiffening coefficient d s . For small extensional strain in the linear range (0 < λ < λ s ), fibers exhibit a constant stiffness k. For compressive strain (− 1 < λ < 0), fibers buckle and show an exponential decay of the stiffness with a characteristic buckling coefficient d 0 . For larger extensional strain beyond the linear range ( λ s < λ ), the fibers show strain stiffening with an exponential increase of stiffness with a characteristic stiffening coefficient d s . Dashed lines indicate the different regimes. Figure adapted with permission from ref. 12 , Springer America Inc.

Extended Data Fig. 2 Rheology and material model for collagen Batch A, B, C.

Collagen hydrogels from three different batches ( a , b , c ) at different concentrations of 0.6 mg/ml (purple), 1.2 mg/ml (blue), and 2.4 mg/ml (orange) are characterized using cone-plate shear-rheological measurements (top; inset shows schematic of experimental setup) and uniaxial stretch-experiments (bottom; inset shows schematic of experimental setup) as described in the Methods section of the main text. The solid lines represent the mean value and the shaded areas represent the standard error for N individual samples. Dashed lines show the fit of the finite element model to the averaged data per condition. The red axis label in c indicates different scaling compared to a , b . d , Material model parameters (Extended Data Fig. 1 ) for non-linear elastic biopolymers (linear stiffness k 0 , the buckling coefficient d 0 , the characteristic strain λ s , and the stiffening coefficient d s .) for each condition. Previously measured parameters from Steinwachs et al. 12 are given for comparison. For each batch, we performed both individual fits to the data for each concentration separately, or global fits where the parameters λ s , d s , and d 0 were the same for all or some of the concentrations. Global fit parameters were preferred if the fit quality was comparable to individual fits, in order to reduce the number of free fit parameters. The increased stiffness of collagen batch C compared to batch A and B is consistent with the observation of decreased cell-generated deformations (Extended Data Fig. 5 ). These stiffness differences can be attributed to batch-to-batch variation. Model parameters for 1.2 mg/ml collagen gels of Batch A were previously published in ref. 58 .

Source data .

Extended Data Fig. 3 Small amplitude rheology.

a , Storage modulus derived from frequency sweeps with a cone-plate shear-rheometer at 1% strain amplitude for different collagen concentrations (Batch A, see Extended Data Fig. 2 ). Dashed lines indicate mean values, and shaded area indicate } one standard deviation. N indicates the number of collagen gels, where each gel was measured in an independent experiment. b , Storage modulus (mean value at 0.02 Hz) scales with collagen concentration according to a power-law with exponent of 1.91 in agreement to previously predicted and measured values 69 , 70 . Bars indicate mean }se and dashed line indicates powerlaw fit curve. Individual collagen gels are shown as dots. The sample size is the same as in a. c , Loss tangent δ (loss modulus G ″ divided by storage modulus \({G}^{{\prime} }\) , averaged between 0.02-2 Hz from the data presented in a) remains below 0.2 for all collagen concentrations,indicating predominantly elastic behavior. d , The FE model parameter k 0 (indicating the linear stiffness of the collagen fibers) for different collagen concentrations (see Extended Data Fig. 2 ) increases approximately linearly with the storage modulus \({G}^{{\prime} }\) of the collagen gels (measured at 0.02 Hz at a strain amplitude of 1% as shown in a). The gray dashed line indicates the predication from continuum mechanics, where k 0 =6E, with Young’s modulus E = 2G(1 + v ) and Poisson ratio ν = 0.25 for linear elastic, isotropic fiber networks 12 . Hence, k 0 = 15 G .

Extended Data Fig. 4 Microstructure of collagen networks.

The collagen fiber structure of two different collagen batches ( a - c : Batch A; e - g : Batch C) is imaged using confocal reflection microscopy for three different collagen concentrations (imaged volume of 160x160x200 μ m with voxel-sizes of 0.314x0.314x0.642 μ m). Grayscale images show a single slice of the imaged volume.3D pore diameters are computed from the covering radius transform as described in 71 , 72 . The mean pore diameters for each collagen concentration ( d , h ) are calculated from 8 different regions within an stack (80x80x100 μ m with 0.314x0.314x0.642 μ m voxel-size). The error bars represent the standard deviation of the mean pore diameter between different regions (shown as individual points) of the imaged stack. Mean value and standard deviation of the pore diameters are listed in the table. See SI Video 1 for 3D representations of the collagen fiber networks.

Extended Data Fig. 5 Dependence of immune cell migration and force generation on matrix stiffness and pore size.

Cell contractility ( a , e ), matrix deformations ( b , f ), cell speed ( c , g ), and cell travelled distance ( d , h ) are measured in collagen gels (0.6 mg/ml, 1.2 mg/ml, and 2.4 mg/ml from different collagen batches) with different stiffnesses ( a - d ) and pore sizes ( e - h ). The storage modulus of different collagen gels is measured with a cone-plate shear-rheometer at 0.02 Hz (1% strain amplitude, see Extended Data Fig. 3 ), and the pore size is derived from confocal reflection images(see Extended Data Fig. 4 ). Colored bars and error bars indicate mean }se for n individual cells (black points) from three (pink and orange bars) or four (green and blue bars) independent experiments. * indicates p<0.05 and ** indicates p<0.01 for two-sided t-test with Bonferroni correction 51 . For clarity, the legend, data points, statistical tests, and cell numbers are only shown in the top row. ( a , e ), Maximum contractility of each cell during a 23 min measurement period. ( b , f ), Maximum of the absolute matrix deformation vector (99% percentile) of each cell during a 23 min measurement period. ( c , g ), Mean cell speed during a 23 min measurement period. ( d , h ), Migration distance of cells after 23 min in 1.2 mg/ml collagen gels. Distance is calculated as the diagonal of the smallest rectangle containing the cell trajectory. Here, only trajectories containing at least 20 data points (corresponding to trajectories of at least 19 min duration) are included in the analysis. Cell contractility monotonically increases, and matrix deformation monotonically decreases with matrix stiffness. Migration speed and travelled distance show a maximum response at intermediate pore sizes.

Extended Data Fig. 6 Force reconstitution for different cell types.

For different cell types, cell-generated matrix deformations are measured and then interpolated onto a finite-element mesh (left). The reconstructed matrix deformations (center) and forces (right) are obtained using Saenopy. The force epicenter is shown in pink, and the size of the image stack is indicated in the lower right. a , Human natural killer cell (NK92) embedded in a 1.2 mg/ml collagen gel (Batch A) during a contractile phase (Fig. 2 , SI Videos 6 – 8 ). b , Hepatic stellate cell (human liver fibroblast) in a 1.2 mg/ml collagen gel (Batch C) after 2 days of culture (SI Video 11 ). c , Axon growth cones of a frog retinal ganglion cell embedded in a 1.0 mg/ml collagen gel (Batch D) during a contractile phase (SI Video 10 ). d , Mouse intestinal organoid in a 1.2mg/ml collagen gel (Batch C) after 24 hours. Time-lapse images of organoid contraction and drug-induced relaxation are shown in SI Video 12 .

Extended Data Fig. 7 Saenopy resolves 3D force fields from bright-field image stacks.

Matrix deformations and forces around a NK92 cell embedded in a 1.2 mg/ml collagen gel (Batch C). Bright-field images are acquired with an ASI RAMM microscope (Applied Scientific Instrumentation, Eugene), CMOS-camera (acA4096-30um, Basler, Ahrensburg),and 20x objective (0.7NA, air, Olympus, Tokyo). The cells are kept at 37 ∘ C and 5% CO2 in a stage incubator (Tokai HIT, Fujinomiya). Matrix deformations (left) are calculated from bright-field image stacks (120 × 120 × 120  μ m, with a voxel-size of 0.15 × 0.15 × 2  μ m, dt = 30 s) during a contractile phase (SI Video 13 , at t = 2 min). Cell forces (right) and corresponding matrix deformations (middle) are reconstructed using a regularization parameter of 10 11 . Maximum intensity projected image stacks are shown below the 3D cubes.

Supplementary information

Supplementary information.

Supplementary Figs. 1–14.

Reporting Summary

Supplementary video 1.

Video showing collagen fibre networks of different concentrations (3D representation): 3D representation of collagen fibre networks (batch C) of different concentrations (0.6, 1.2 and 2.4 mg ml −1 ). Collagen fibres were imaged using confocal reflection microscopy. Images are Sato-filtered for ridge detection, highlighting the fibre structure 73 . The recorded stack size is 160 × 160 × 50 μm with a voxel size of 0.314 × 0.314 × 0.642 μm.

Supplementary Video 2

Video of an NK92 cell migrating through collagen (bright-field and confocal reflection): intensity projected image stacks of an NK92 cell migrating through a 1.2 mg ml −1 collagen gel (stack size of 123 × 123 × 123 μm with a voxel size of 0.24 × 0.24 × 1 μm) recorded every minute. The cell shows a phase of high contractility at around t  = 4 min. Left: bright-field image. Right: confocal reflection image. Projected matrix deformations are indicated by coloured arrows. Scale bar, 20 μm. For better visualization of matrix deformations, the video is followed by a forward–backward sequence of consecutive images during a phase with high contractility.

Supplementary Video 3

Video of an NK92 cell migrating through collagen (differential interference contrast): differential contrast images of an NK92 cell migrating through a 1.2 mg ml −1 collagen gel. The time between consecutive frames is 2 seconds. The focus has been adjusted manually during the recording. For better visualization of matrix deformations, the video is followed by a forward–backward sequence of consecutive images during a phase with high contractility.

Supplementary Video 4

Video of an NK92 cell migrating through collagen (3D representation): 3D representation of an NK92 cell (green) migrating through dense constrictions in a 1.2 mg ml −1 collagen gel (brown). The cell was imaged using calcein staining (2 μM calcein AM, Thermo Fisher Scientific) and segmented using Yen thresholding 74 . Collagen fibres were imaged using confocal reflection microscopy. Images are Sato-filtered for ridge detection, highlighting the fibre structure 73 . The transparency of the cell and collagen image stacks were determined by the intensity values according to a sigmoidal transfer function. The recorded stack size is 123 × 123 × 40 μm with a voxel size of 0.24 × 0.24 × 1.5 μm. The time between consecutive image stacks is 1 min. Matrix deformations are indicated by coloured arrows. For better visualization of matrix deformations, the video is followed by a forward–backward sequence of consecutive images during a phase with high contractility.

Supplementary Video 5

Video of an NK92 cell migrating through collagen (3D representation): 3D representation of an NK92 cell (green) migrating through dense constrictions in a 1.2 mg ml −1 collagen gel (brown). The cell was imaged using calcein staining (2 μM calcein AM, Thermo Fisher Scientific) and segmented using Yen thresholding 74 . Collagen fibres were imaged using confocal reflection microscopy. Images are Sato-filtered for ridge detection, highlighting the fibre structure 73 . The transparency of the cell and collagen image stacks were determined by the intensity values according to a sigmoidal transfer function. The recorded stack size is 123 × 123 × 40 μm with a voxel size of 0.24 × 0.24 × 1.5 μm. The time between consecutive image stacks is 1 min. Matrix deformations are indicated by coloured arrows.

Supplementary Video 6

Video of matrix deformations around an NK92 cell migrating through collagen: deformation field around an NK92 cell during migration in a 1.2 mg ml −1 collagen gel (stack size of 123 × 123 × 123 μm with a voxel size of 0.24 × 0.24 × 1 μm, d t  = 1 min). The cell shows a phase of high contractility at around t  = 5 min. Intensity projected bright-field image stacks are shown below the 3D fields. Note that the contractile phase occurs when the cell moves through a narrow constriction.

Supplementary Video 7

Video of matrix deformations calculated from confocal reflection images of collagen fibres around an NK92 cell: the deformation field around the NK92 cell shown in Supplementary Video 6 at t  = 5 min (peak of the contractile phase). The deformation field is calculated from confocal reflection images using PIV 27 . Individual images are recorded using a resonance scanner operated at 8,000 Hz in combination with a galvo stage to move the sample in z direction. The entire image stack is recorded within 10 seconds (no line or frame averaging). Despite substantial image noise, the deformation field can be reliably calculated.

Supplementary Video 8

Video showing the force reconstruction of a migrating NK92 cell: measured matrix deformations (left), reconstructed matrix deformations (centre) and reconstructed force field (right) around the NK92 cell shown in Supplementary Video 6 at t  = 5 min (peak of the contractile phase). The measured matrix deformations agree well with the reconstructed deformations. Intensity projected bright-field image stacks are shown below the 3D fields. The pink dot denotes the force epicentre of the force field.

Supplementary Video 9

Video showing elastic matrix deformations around NK92 cells: NK92 cells embedded in 1.2 mg ml −1 collagen gels during contractile phases imaged by confocal reflection microscopy. To visualize changes in collagen structure, the initial state of the matrix is shown in red and the current timestep is shown in green. The video is followed by a sequence of images showing the state of the matrix before, during and after a contractile burst. The collagen fibres relax back to their initial position, indicating predominantly elastic matrix deformations without notable plastic deformations.

Supplementary Video 10

Video of axon growth cones: maximum projected confocal reflection image stacks of Xenopus retinal ganglion cell axon growth cones in a 1.0 mg ml −1 collagen gel (batch D). The recorded stack size is 145 × 145 × 50 μm with a voxel size of 0.14 × 0.14 × 1 μm. The time between consecutive image stacks is t  = 5 min. Similar to contractile phases of immune cells, we observe contractile phases during axon growth cone development. For better visualization, we show a projected volume of 124 × 124 × 10 μm around the axon growth cones and show a forward–backward sequence of the contractile phases at the end of the video.

Supplementary Video 11

Videos showing force reconstruction of a hepatic stellate cell: confocal reflection and fluorescence image stacks around a human hepatic stellate cell in a 1.2 mg ml −1 collagen gels are acquired after a culture time of 2 days (stack size of (370 μm) 3 with a voxel size of 0.72 × 0.72 × 0.99 μm). Matrix deformations were calculated from the confocal reflection image stacks before and after relaxation of cellular forces using 10 μM cytochalsin D. The measured matrix deformation field (left) agrees well with the reconstructed deformation field (centre). Intensity projected image stacks of the calcein-stained cell are shown below the 3D fields. The purple outline indicates the 3D reconstruction of the calcein-stained cell and the pink dot represents the force epicentre (right).

Supplementary Video 12

Video of an intestinal organoid in collagen: image stacks of an intestinal organoid in collagen (1.2 mg ml −1 ) were recorded using confocal reflection microscopy (stack size 738 × 738 × 100 μm with a voxel size of 0.74 × 0.74 × 2 μm, d t  = 20 min). Organoids were relaxed after 24 hours using 10 μM cytochalsin D and 0.1% Triton x-100. Images are Sato-ridge filtered to highlight the fibre structure 73 . The full time-lapse sequence is followed by a forward–backward sequence of an image pair taken before and after relaxation of cellular forces to visualize the total contraction.

Supplementary Video 13

Video showing a contractile phase of a migrating NK92 cell measured from bright-field image stacks: bright-field image stacks of a NK92 cell during migration in 1.2 mg ml −1 collagen gel (batch C) were acquired with an ASI RAMM microscope (Applied Scientific Instrumentation), a CMOS camera (acA4096-30um, Basler) and a ×20 objective (0.7 numerical aperture, air, Olympus). Cells were kept at 37 °C and 5% CO 2 in a stage incubator (Tokai HIT). Matrix deformations were calculated from bright-field image stacks (d t  = 30 s, volume 120 × 120 × 120 μm, with a voxel size of 0.15 × 0.15 × 2 μm). Contractions are visualized by a forward–backward sequence of an image pair taken 30 s apart. The left panel shows the maximum intensity projected bright-field image stack (40 μm height around the centred cell) with differential matrix deformations (arrows). The matrix deformations are clearly visible even though the individual collagen fibres cannot be resolved (right, single image plane).

Source data

Source data fig. 1.

Source data.

Source Data Fig. 2

Source data fig. 3, source data fig. 4, source data fig. 5, source data extended data fig. 2, source data extended data fig. 3, source data extended data fig. 4, source data extended data fig. 5, rights and permissions.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions

About this article

Cite this article.

Böhringer, D., Cóndor, M., Bischof, L. et al. Dynamic traction force measurements of migrating immune cells in 3D biopolymer matrices. Nat. Phys. (2024). https://doi.org/10.1038/s41567-024-02632-8

Download citation

Received : 08 September 2023

Accepted : 06 August 2024

Published : 24 September 2024

DOI : https://doi.org/10.1038/s41567-024-02632-8

Share this article

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

Quick links

• Explore articles by subject
• Guide to authors
• Editorial policies

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.