High-speed photoelastic tomography for axisymmetric stress fields in a soft material: temporal evolution of all stress components

This study presents a novel approach for reconstructing all stress components of the dynamic axisymmetric fields of a soft material using photoelastic tomography (PT) and a high-speed polarization camera. This study focuses on static and dynamic Hertzian contact as an example of transient stress field reconstructions. For the static Hertzian contact (a solid sphere pressed against a gel block), all stress components in the urethane gel, which has an elastic modulus of 47.4 kPa, were reconstructed by PT using the measured photoelastic parameters. The results were compared with theoretical solutions and showed good agreement. For the dynamic Hertzian contact (a sphere impacting gel), a high-speed polarization camera was used to reconstruct the transient stress field within the gel. PT was used to quantitatively measure the shear and axial stress waves and showed different propagation speeds on the substrate. The technique allowed the simultaneous measurement of stress fields ranging from $O(10^{-1})$ to $O(10^1)$ kPa during large deformations, demonstrating its accuracy in capturing rapidly changing stress tensor components in dynamic scenarios. The scaling laws of the calculated impact force agreed with theoretical predictions, validating the accuracy of PT for measuring dynamic axisymmetric stress fields in soft materials.


Introduction
Evaluating stress in soft materials is of great importance in the fields of biomedical engineering and cell printing applications.For example, in biomedical engineering, it is essential to determine the stress distributions around cerebral aneurysms to decipher the mechanisms leading to their rupture [1][2][3].Investigating the stress field induced by needle-free injection for evaluating pain levels also requires stress measurements in the tissue, which is also a soft material [4,5].Moreover, understanding the stress distribution within a material resulting from droplet or liquid jet impact is pertinent to diverse engineering processes, making it a topic of considerable interest for research [6][7][8][9][10].
Photoelasticity is a well-known technique for estimating stress distributions within materials [11][12][13][14] and is expected to be useful in the scenarios mentioned above.Photoelasticity operates on the principle of birefringence in materials.Birefringence is an optical anisotropy, in which the local refractive index of a stressed material changes in different directions.By observing the intensity of polarized light passing through the stressed body, one can measure the difference between two mutually orthogonal refractive indices and the orientation of its axis to estimate the material's internal stress state.The stress-optic law establishes a proportional relationship between the integration of the secondary principal stress difference σ d in a material and the phase retardation of transmitted light ∆ [12,15], ∆ = C σ d (y)dy.Here, the coefficient C -known as the stress-optic coefficient -is specific to the material [12], and y is the camera's optical axis.The secondary principal stress is the apparent principal stress that appears in the cross-section perpendicular to y-axis [16].Additionally, the orientation ϕ of the transmitted polarized light is related to the direction of the secondary principal stress.In this study, the retardation and orientation of the polarized light are called photoelastic parameters.
Photoelasticity possesses distinct advantages over other measurement methods, such as pressure sensors, due to its ability to measure the entire stress field within a material while being non-invasive.Digital image correlation (DIC) [7,17] can also be used to obtain a full-field stress measurement by analyzing local material displacements and applying the material's constitutive equation.Such a technique requires particles to be mixed into the material as a marker to measure the displacement or optics and to produce very thin laser sheets.In contrast, photoelasticity can obtain measurement data on stress values without contaminating the material by installing markers, as long as the stress-optic law is valid.Additionally, the experimental setup of photoelasticity requires only circularly polarized light as a light source and a polarization camera.
Methods designed for measuring three-dimensional stress fields in materials are called integrated photoelasticity methods [14,[18][19][20].In this context, "integrated" signifies that the recorded retardation captures incremental photoelastic parameters at each material point along the path of the light ray, depending on the stress state at that specific point.Previously, integrated photoelasticity has mainly been used for hard materials such as glass to visualize its residual stress [11-14, 21, 22].In contrast, our previous study verified the applicability of integrated photoelasticity to a soft material experiencing an axisymmetric stress field [16,23].Consequently, we have demonstrated that integrated photoelasticity is valid for a gel that is about 1,000,000 times softer than glass.
As mentioned above, the photoelastic parameters (the phase retardation and orientation) of polarized light passing through a stressed material are related to the integrated values of the stress field inside the material, as shown in Fig. 1.Therefore, it is necessary to reconstruct the internal stress field from the measured photoelastic parameters to obtain the stress field in a material.This is called photoelastic tomography and was studied intensively by Aben et al. [24][25][26].Tomography is well known as a nondestructive method for analyzing the internal structure of three-dimensional objects.For example, computed tomography (CT) scans, which are based on X-rays, use the attenuation of the X-rays transmitted through an object to determine its internal structure.In such cases, tomography is relatively easy because the quantities that characterize the internal structure at each point are scalar values.On the other hand, stress is a tensor with up to six components at each point inside the material.Two values are obtained by integrating them -the phase retardation and orientation -which is a small number compared to the number of tensor components to be reconstructed.Therefore, tomography of such a field is generally difficult, but if the stress field is axisymmetric, tomography of the stress field can be performed because the number of tensor components is reduced [27,28], i.e., Here, σ rr , σ θθ , σ zz , and σ rz are the radial, hoop, axial, and shear stresses in cylindrical coordinates, respectively.However, photoelastic tomography has not yet been applied to the full-field measurement of the internal stress of a soft material.Additionally, photoelasticity has been used primarily to measure steady-state stress fields, because measuring photoelastic parameters using the phase-shifting method requires multiple images to be taken for the same stress field.However, the development of a polarization camera with polarizers of different angles installed in the camera sensor has made it possible to measure photoelastic parameters in a single shot [29].Furthermore, highspeed polarization [30] can measure transient photoelastic parameters even for very high-speed phenomena.However, although the high-speed polarization camera has been used to "visualize" dynamic three-dimensional stress fields (i.e., to measure photoelastic parameters) in a soft material [31], to the best of our knowledge, there have been no examples of quantitative measurement through the reconstruction of stresses in soft materials from the measured photoelastic parameters based on photoelastic tomography.
In the present work, photoelastic tomography is applied to measure all stress components of dynamic axisymmetric fields in a soft material using a high-speed polarization camera to address the gaps described above.Combining photoelastic tomography and the high-speed polarization camera allows us to reconstruct all the components of the dynamic axisymmetric stress fields in a soft material.The ability to measure all components of the stress tensor not only allows analysis of the principal stresses and their directions (rather than the secondary principal stresses) but also allows the measurement of the von Mises stress, which is used to evaluate the yield points of the material [23].The examination conducted in this study involves both static and dynamic Hertzian contact problems.There are two primary reasons for this.Firstly, if the external force and other contact conditions are known a priori, the stress field within the material or impact forces acting on the surface of the material can be obtained theoretically [32][33][34][35][36].This allows us to evaluate the certainty of the reconstructed stress field from the experimentally obtained photoelastic parameters.Secondly, dynamic axisymmetric stress fields appear in many phenomena, e.g., a droplet impacting a soft surface [37,38], a liquid jet penetrating a gel [4,5,31], and cavitation generation in a gel [39].For these two reasons, this study demonstrates that photoelastic tomography is an effective technique for measuring dynamic axisymmetric stress fields in materials.We believe that photoelastic tomography will find broad applications as a method for understanding the high-speed phenomena of soft materials.
In Sec. 2, the integrated photoelasticity and the photoelastic tomography algorithm are described.The experimental method used to measure the photoelastic parameters with a high-speed polarization camera is also explained.Section 3 compares the experimental and theoretical results, and Sec. 4 summarizes the findings of this study.

Integrated photoelasticity
When circularly polarized light from a light source enters a material under stress, it is emitted as elliptically polarized light with photoelastic parameters (phase retardation ∆ and orientation ϕ) corresponding to the stress state in the material [18].∆ and ϕ can be temporally measured using a high-speed polarization camera [30].The following relationship exists between the stress field inside a material and its photoelastic parameters [40][41][42]: where C is a material-specific value called the stress-optic coefficient and σ xx , σ zz , and σ xz are the stress components in Cartesian coordinates with the y-axis as the camera's optical axis.If ∆ is smaller than λ/4 and the rotation of the principal stress direction does not exceed π/6, Eqs. ( 2) and (3) are valid [12,19].Note that all of the experimental data in this study are within the bounds of this assumption.As shown by the integrated equations, the photoelastic parameters that can be measured using a polarization camera are the integrated values of the internal stress field along the camera's optical axis.Therefore, tomography is required to convert the photoelastic parameters into the stress tensor components at each point inside the material.However, as mentioned in the previous section, the tomography of stress fields differs from traditional tomography, which targets a problem where a single scalar characterizes each point inside the material.Because stress is represented by a tensor, applying conventional scalar tomography directly to a stress field is difficult.Aben et al. suggested that the problem of stress field tomography could be solved if it could be reduced to a problem of scalar field tomography for a single stress tensor component [18,24].To determine the stress field inside a material using photoelastic tomography, it is necessary to measure the photoelastic parameters in two sections, the upper section and lower section, parallel to the x-y plane and separated by a distance ∆z [28,43].Generally, this measurement should be implemented for many directions of the camera's optical axis around the z-axis.We now consider an arbitrary three-dimensional stress field with boundaries (Fig. 2(a)).In a part of this field, the element ABX (see Fig. 2(b)), the force equilibrium along the x-axis is ∆z where T u and T l are shear forces on the upper and lower surfaces of the element and A, B denote the boundary of the three-dimensional stress field on the y-axis.T u and T l can be described using Eq.(3): Here, V ′ 2 is the V 2 measured in the upper section, and X ′ and X are the outer positions of the upper and lower sections, respectively.Substituting Eqs. ( 4) and (5) into Eqs.
(2) and (3), we obtain Using this equation, the tensor tomography of the stress field can be treated as scalar field tomography.Note that the reconstruction of the axial stress using Eq. ( 7) requires the assumption of the equilibrium of the force, as mentioned above.In contrast, the reconstruction of the shear stress using Eq. ( 6) does not require this assumption.

Photoelastic tomography
If the stress field is axisymmetric, for example, in the problem of a rigid sphere making contact perpendicular to a flat surface, the stress field can be experimentally determined using Eqs.( 6) and ( 7) [27,28].When the stress field is axisymmetric, the stress tensor has four components: radial σ rr , hoop σ θθ , axial σ zz , and shear stresses σ rz , i.e., Eq. ( 1).Here, the plane perpendicular to the camera's optical axis is the r-z plane, and the axis of symmetry is the z-axis.The components σ rz and σ zz are reconstructed from the measured values of ∆ and ϕ using the integrated relationships provided in the last section.The remaining components, σ rr and σ θθ , are reconstructed from σ rz and σ zz using the equations of the theory of linear elasticity [27].In this subsection, tomography for general axisymmetric scalar fields and its algorithm will be described, followed by the application of the algorithm to tomography for stress fields.

Abel transform and the onion-peeling method
A well-known reconstruction method for an axisymmetric field is the inverse Abel transform.Generally, the relationship between an axisymmetric field p(r) and a projection function P (x), as shown in Fig. 3, which is obtained by integrating over p(r), is called Abel transform and can be expressed as: where r = x 2 + y 2 .This equation can be rewritten as The inverse Abel transform is given by Since the projection data obtained in the experiment are always discretized, this inverse transformation must be solved numerically.Direct use of the above equation not only gives rise to a singularity at x = r, but also requires that the exact change in P in the x direction be known a priori.However, the projection data obtained in the experiment are not perfectly smooth.It should also be noted that Abel's forward and inverse transformations cover the entire space and have infinite integration limits, whereas the experimental data are limited to a finite range.Outside of these ranges, p and P must be strictly zero, and any deviation will compromise the accuracy of the transform.Basically, only local objects that are completely contained within the image frame can be accurately transformed.The onion-peeling (OP) method is one of the numerical algorithms used to solve the inverse transform [44,45].The OP method considers the axisymmetric field as layers of N rings (Fig. 3), which have a thickness of ∆r and a fixed scalar value.The outermost ring's value p(r N ) is obtained from the projection data P (r N ) and the length of the optical path 2W N,N using the equation The value of the N -1th ring p(r N −1 ) can be obtained by using the projection data P (r N −1 ), the value of the N th ring p(r N −1 ) obtained earlier, and the length of the optical path for each ring, 2W N −1,N −1 and 2W N −1,N , i.e., Therefore, the projection data of the ith ring can be described as where r i = i∆r is the distance from the central axis and

Determination of the shear and axial stresses
We determine the shear stress σ rz from the photoelastic parameters using Eq. ( 3) based on the OP method.While the OP method described earlier targets axisymmetric fields, the shear stress σ rz is no longer axisymmetric when it is distributed in the x-y plane in Cartesian coordinates.The relationship between the shear stresses in Cartesian coordinates σ xz and cylindrical coordinates σ rz is Therefore, the angle θ has to be considered when we relate the integrated value of the photoelastic parameters and the shear stress.This is the main difference from the simple onion-peeling method described previously.The stress σ rz to be reconstructed and the stress σ xz to be integrated are linked by an angle θ, as in Eq. ( 15) [28].This is because photoelasticity can measure only the stress component projected on the plane perpendicular to the camera's optical axis.In the N th ring, is the integral of the shear stress σ (N ) rz acting on the N th ring, i.e., where W N,N is the light path of the N th camera's optical axis passing through the N th ring.From this equation, the shear stress σ rz , can be determined from the following equation: In the same way, the shear stress of the ith ring can be expressed as where The axial stress σ zz can be determined similarly.From Eq. ( 7), the axial stress acting on the ith ring σ zz can be determined using the following equation: In the literature [28], it has been shown that the stresses can be calculated recursively from the outermost ring inward using Eqs.( 18) and ( 20), and finally, the shear and axial stresses for all rings can be calculated.
Here, we transform the summation equations, Eqs. ( 18) and ( 20), into a matrix form represented by the product of the stress component vector, the coefficient matrix, and the integrated value vector [46].Equation ( 18) becomes where (2) 2 . . .
The coefficient matrix for the shear stress α rz is described as where α i,j = W i,j cos θ i,j .
In the same way, the summation equation for the axial stress, Eq. ( 20), can be rewritten as where . . .
The vector of the shear and axial stresses σ rz , σ zz can be obtained instantaneously by multiplying the vector of the integrated value by the inverse of the coefficient matrices, i.e., if the inverse matrices, α −1 rz and α −1 zz , are calculated in advance.This reduces the computational time required to calculate stresses acting on all rings compared with the recursive calculation of Eqs.(18) and (20) from the outermost ring inward.

Determination of the radial and hoop stresses
The radial σ rr and hoop σ θθ stresses are determined using the equations of the theory of linear elasticity [25,27].The stresses must satisfy the equilibrium equation and compatibility equation, respectively.The differential equations for the radial and hoop stresses are obtained by transforming these equations as follows: Since the spatial distribution of the shear and axial stresses was obtained in Sec.2.2.2, the second and subsequent terms on the right-hand side of these equations are known values.Therefore, by solving this differential equation numerically, the radial and hoop stresses can be determined.Dis-cretizing Eqs. ( 32) and ( 33) using the first-order forwarddifference method yields θθ =σ −ν ∂σ rz ∂z where i = 1 for r = 0 and i increases toward r.The radial and hoop stresses can be determined by setting appropriate boundary conditions.

Quasi-static approximation
The tomography algorithm described in the previous sections assumed equilibrium except for the case of the reconstruction of the shear stress σ rz .Therefore, the quasi-static approximation has to be adopted to reconstruct stresses from time-series data of the measured photoelastic parameters.The impact velocity V and the pressure wave speed v p in the material are compared to confirm this approximation.The pressure wave speed v p can be expressed using the following equation [32]: where E and ρ s are Young's modulus and the density of the substrate material, respectively.For low-velocity impacts, where the sphere velocity V is much smaller than the pressure wave velocity v p , the quasi-static approximation is more likely to be valid.

Experiment
The schematics of the setups for the experimental verification in the static and dynamic cases are shown in Figs.4(a) and 4(b), respectively.The setup includes a light source, a measurement target, and a polarization camera.A urethane gel block (urethane gel phantom, Exseal Co., Ltd., Gifu, Japan), which has a size of 50 × 50 × 50 mm 3 and a density of ρ s of 1064 kg/m 3 , is placed on an electronic balance.We estimated the Young's modulus of the gel, E, from the surface deformation δz max along the z-axis and the loading normal force F n with the electric balance using Hertzian contact theory [32], i.e., where E * = E/(1 − ν 2 ).The Young's modulus of the gel was determined to be 47.4 kPa, with the Poisson's ratio of the gel ν assumed to be 0.5 in this study.Note that Hertzian contact theory assumes that the contact radius (the area where the sphere touches the substrate) is considerably smaller than R [32,47].However, the measured surface deformation results closely follow the curve outlined in Eq. ( 37) (see Fig. 5).Our previous paper delves into the Hertzian contact problem in the context of a highly deformable substrate [48].
For the static case (Fig. 4(a)), a styrol sphere with a diameter of 2R = 14.7 mm is vertically pressed against the surface of the gelatin with a loading normal force F n .This experimental setup produces Hertzian contact between a sphere and a half-space [32].The loading normal force F n is calculated as F n = mg, where the apparent mass m is measured using the electronic balance.The acceleration of gravity used in the calculation is 9.81 m/s 2 .The retardation and orientation fields under sphere pressing are recorded by the polarization camera (Photron, CRYSTA PI-5WP) with the light source (Thorlabs, SOLIS-565C).To reduce the deterioration of the measurement accuracy due to the bandwidth of the light source wavelength [49], a narrow bandwidth bandpass filter (Photron, with a half-bandwidth of approximately 20 nm) is installed in front of the camera sensors.Then, the typical wavelength of the light is 540 nm.The number of pixels in the measured images of the photoelastic parameters is 191 × 208 pixels in the region of 1.37R × 1.50R, which is used for stress reconstruction.The spatial resolution of the photoelastic parameters images is 52.7 µm/pix for the static case.
For the dynamic case (Fig. 4(b)), a plastic sphere with a radius of R = 2.98 mm and a density of ρ = 997.3kg/m 3 impacts the gel during free fall.The impact velocity of the sphere, V , is varied from 0.3 to 0.7 m/s by adjusting the falling height of the sphere.The retardation and orientation fields during the sphere impact are recorded by the polarization camera (Photron, CRYSTA PI-1P) with the light source (Thorlabs, SOLIS-525C).The typical wavelength of the light is set to 520 nm by using a bandpass filter placed in front of the lens attached to the camera.The temporal resolution of the camera is 20,000 f.p.s for high-speed stress measurement.The number of pixels in the measured images of the photoelastic parameters is 115×171 pixels in the region of 3.28R × 5.14R, which is used for stress reconstruction.The spatial resolution of the photoelastic parameter images is 84.9 µm/pix for the dynamic case.
The light source generates incident light of a typical wavelength λ.After passing through the polarizer, which has an angle of 0 • , the incident light then goes through a quarterwave plate, which has an angle of 45 • , resulting in circular polarization.The stressed gel allows circularly polarized light to pass through and emits it as elliptically polarized light with retardation ∆ and orientation ϕ.Using a polarization camera, it is possible to simultaneously measure both the orientation and retardation of the emitted light.Four linear polarizers with angles of I 0 • , I 45 • , I 90 • , and I 135 • are installed in the four neighboring pixels of the image sensor of the polarization camera.This set of four neighboring pixels is called a "super-pixel" [50].The intensity values measured by the camera's sensor through each linear polar-  Using the four-step phase-shifting method [15,16,30,51], the retardation ∆ and orientation ϕ are obtained from the four intensity values of the super-pixel as follows: where As an example of the measured data, the retardation and orientation are calculated using software (Photron Ltd., CRYSTA Stress Viewer) and are shown in Figs.6(a) and 6(b).The retardation and orientation data are filtered to reduce the noise using Matlab's median filter function (medfilt2).The details of the principle of measurement are also described in Refs.[16,30,52].
The theoretical retardation and orientation field can be obtained using the integrated photoelasticity [16] from the theoretical solution of the stress fields calculated using the Hertzian contact problem [16,23,32,33], for which an example of the theoretical data is shown in Figs.6(c) and 6(d).The retardation and orientation fields are calculated using the same conditions as in the experiment shown in Figs.6(a) and 6(b) using a stress-optic coefficient of 1.14 × 10 −9 Pa −1 .The stress-optic coefficient was deter-mined by comparing the experimental and theoretical retardation fields, and chosen to be the value that minimizes the root mean square error between the calculated and experimental retardation.The procedure for determining C is described in our previous paper [16].The result of the determination of the stress-optic coefficient is shown in Fig. 7.The stress-optic coefficient is almost constant over a wide range of loading normal forces F n .

Experimental validation of photoelastic tomography for static Hertzian contact
In this section, we demonstrate that photoelastic tomography can reconstruct all the components of the axisymmetric stress fields, σ rr , σ θθ , σ zz , and σ rz , in a soft material that shows large deformation, focusing on the static Hertzian contact problem where a rigid sphere is pressed against an elastic half-space.Figure 8 compares the reconstructed stress fields (righthand side), which are from the experimental photoelastic parameter data, and the theoretical stress fields (left-hand side), which are obtained by solving the Hertzian contact problem using similar conditions to the experiment for m, R, and E. For all stress components, there are errors near the surface of the gel, but overall, the experimental and theoretical results agree well.It should be noted that the experimental results do not necessarily completely agree with the theoretical results because the contact conditions between the sphere and the gel surface are not necessarily identical in the theoretical calculations and the experiments.In other words, errors between the experimental and theoretical results are expected to be a result of the surface contact conditions in the experiment, which are not considered in the calculation of Hertzian contact theory.The following issues can be considered as the sources of the noisy distribution of the shear stress σ rz around the z-axis.The orientation values approach zero around the z-axis (see Fig. 6(b)) and are sometimes negative.Theoretically, the azimuth should be greater than or equal to zero on the right-hand side of the z-axis.Therefore, the filtering of the orientation data before reconstruction and the determination of the z-axis position has a significant influence on this noise.Additionally, the axial stress σ zz distribution shows  horizontal stripe noise, while the shear stress σ rz does not show any.This is because the z-axis gradient of the V 2 field, i.e., V ′ 2 − V 2 , is used when reconstructing the axial stress (see Eq. ( 20)).Since the V 2 distribution is not perfectly smooth, the V ′ 2 − V 2 distribution is discontinuous along the z-axis.Therefore, the axial stress calculated using this distribution is also not smooth in the z-direction, resulting in horizontal stripe noise.This noise carries over to the radial and hoop stresses.
Since all stress components have now been determined, the von Mises stress σ M , which is often used as a failure criterion for ductile materials and as a force criterion related to plastic deformation [23], can also be determined using the equation where Figure 8(e) compares the theoretical and experimental von Mises stresses and shows good agreement.
Figure 9 shows the line profiles of the reconstructed stress fields compared with the theoretical solutions at a certain z position, z/R = 0.7.The reconstructed stress fields show similar distributions to those estimated by solving the Hertzian contact problem.The noisy distribution of the line profiles is mainly due to the noise of the photoelastic parameters (the retardation and orientation fields).In particular, in the calculation of the axial stress σ zz , the integral value vector expressed in Eq. ( 27) is sensitive to the noise of the photoelastic parameters because it is the amount of change between adjacent pixels in the z-direction of V 2 consisting of photoelastic parameters.Smoothing the experimental data with a more appropriate filter is expected to reduce this noise.
These results indicate that photoelastic tomography can reconstruct all components of axisymmetric stress fields in soft materials with large deformation.This allows for the estimation of the von Mises stresses, which are important for evaluating the yield value of a material, and the analysis of principal stresses and their directions (rather than the secondary principal stress).

High-speed stress field measurement during dynamic Hertzian contact
In this section, we determine the applicability of photoelastic tomography for high-speed stress field measurements in soft materials, focusing on the dynamic Hertzian contact problem, in which a rigid sphere impacts an elastic half-space.Firstly, the validity of the quasi-static approximation in this dynamic case was checked by comparing the sphere's impact velocity V and the substrate's pressure wave speed v p (Sec.2.2.4).In our experiment, the maximum value of V /v p is about 0.4 because the maximum V is 2.7 m/s while v p is approximately 6.7 m/s, since E = 47.4 kPa and ρ s = 1064 kg/m 3 .Therefore, we assumed that the quasistatic approximation holds in our experiment since the ratio V /v p is smaller than 1.
Figure 10 shows the retardation and orientation distribution in a gel when a plastic sphere with a radius R = 2.9 mm makes impact at a velocity V ≃ 2.2 m/s.The retardation and orientation corresponding to the dynamic stress field in the soft material caused by the sphere impact are measured.The white area indicates the gel's surface de-formation due to the sphere's penetration.Before impact (t = 0 ms), because there is no stress in the gel, the retardation is zero, and the orientation is random from −90 • to 90 • .The maximum retardation can be observed at around t = 0.5 ∼ 0.75 ms, when the vertical displacement of the sphere reaches its maximum value.Then, the retardation decreases with time according to the upward motion of the sphere.
Figures 11 and 12 show the temporal evolution of the shear, axial, radial, and hoop stress distributions during the impact of the sphere, as reconstructed from the photoelastic parameters (Fig. 10).The accuracy of the reconstructed stress fields was validated by comparing the scaling law for the maximum impact force F n,max calculated from the axial stress field with the theoretical prediction.The details of this comparison are discussed in the next section (Sec.3.3).A stress maximum appears between 0.25 and 0.5 ms for all stress components.After that, the stresses decrease with time following the upward motion of the sphere.For the shear and axial stresses, negative values eventually appear near the edges of the contact area, probably due to the adhesion force between the spheres and the gel surface.For  the radial stress, the region of the acting stress is wider than the others, and a negative stress region appears along the z-direction.The hoop stress exhibits both negative and positive values immediately following the impact, and then the negative stress is damped.
Figure 13 shows kymographs of the shear and axial stresses acting on the substrate surface, σ rz (r, z = 0, t) and σ zz (r, z = 0, t).Figures 13(b) and 13(d) are plotted using the same data as in Figs.13(a) and 13(c), respectively, with different color bar limits to visualize stress values with small magnitudes.The waves of the shear and axial stresses at the surface show different velocities.This cannot be obtained from the visualization of the retardation and orientation data measured by the polarization camera alone, such as in Fig. 10.It can only be quantified by reconstructing the stress field using these data.Figures 11,12,and 13 show that even stress fields with a wide range from O(10 −1 ) to O(10 1 ) kPa can be measured simultaneously.The measurable range of stress depends on the wavelength of the light source λ and the stress-optic coefficient of the material used, C. From Eq. ( 38), the measurement range of retardation ∆ is limited to 0 −λ/4.Also, from Eqs. ( 2) and ( 3), the magnitude of ∆ is determined by the stress-optic coefficient C. Therefore, the order of magnitude of the measurable stress can be changed by using a different gel, i.e., changing the stress-optic coefficient C. Gelatin gel has an order of magnitude larger C [12,16,53,54] than the urethane gel used in this experiment and is more sensitive to small stresses.Additionally, the range of measurable retardation could be further expanded by using a color polarization camera [55].
These results demonstrate that photoelastic tomography can be applied successfully even in fields where the stress tensor components are fast and time-varying, such as stress wave propagation.

The scaling law for the maximum impact force
In this section, the normal impact force acting on the gel surface is estimated, and the scaling laws for its maximum values are compared with theoretical predictions from dynamic Hertzian contact to evaluate the accuracy of the photoelastic tomography in reconstructing the time-varying stress fields.
The normal impact force F n can be obtained by integrating over the axial stress σ zz acting on the surface of the gel, i.e., Figure 14 shows the temporal evolution of the normal impact force F n for different impact velocities.The normal impact force suddenly increases after impact and reaches its maximum value.Then, it decreases following the sphere's vertical motion because of the gel's viscous damping effect.The maximum normal impact force increases with the impact velocity.In the case of a high impact velocity, the impact force becomes negative when the sphere moves upward.This is likely due to the adhesion force between spheres and the gel surface.
Here, we focus on the peak force F n,max and the time when the peak force occurs t Fn,max .The scaling law of the peak force and time with the impact velocity V can be theoretically predicted from the dynamic Hertzian contact model [35,36].The motion of the equation of a sphere impacting the elastic half-space can be described as with initial conditions dz/dt = V and z = 0 at t = 0, where z is the vertical displacement of the sphere.Here, we neglect the gravity term because it is negligibly small compared with the other terms.The scaling law of F n can be theoretically obtained from the Hertzian contact problem using the displacement of the surface [32], which corresponds to z as where E * = E/(1 − ν 2 ) is the effective Young's modulus of the gel.Here, by defining v = dz/dt, Eq. ( 43) can be rewritten using Eq. ( 44) as Thus, Eq. ( 45) can be solved for v in the following equation.
Because v = 0 when z reaches its maximum z max , By substituting this equation into Eq.( 44), the scaling law for the peak force F n,max can be written as The peak force F n,max can be nondimensionalized using the inertial force, ρV 2 R 2 , as   Additionally, we can obtain the scaling law for the peak time t Fn,max from the relationship z max ∝ t Fn,max V and Eq. ( 47), as t Fn,max ∝ E −2/5 ρ 2/5 RV −1/5 . ( Let this equation be nondimensionalized by the impact Figure 15(b) shows the relationship between t Fn,max /T and ρV 2 /E with Eq. ( 51) as a dashed line.The trend of the experimental peak time agrees with the theoretical prediction.
Note that Hertzian contact theory, including the dynamic cases, basically assumes a small deformation of the substrate surface.Nevertheless, the theoretical predictions of the scaling law still agree well with the experimental data over a wide range of ρV 2 /E from O(10 −3 ) to O(10 −1 ), even with large substrate deformations.

Conclusion
In this study, all stress components in the dynamic axisymmetric fields of a soft material have been reconstructed using photoelastic tomography.We have validated our approach, focusing on static and dynamic Hertzian contact.
For static Hertzian contact, in which a rigid sphere is pressed against a gel block, all stress tensor components were reconstructed from the photoelastic parameters measured using a polarization camera.The reconstructed stress fields were compared with the theoretical solution for the stress fields obtained by solving the Hertzian contact problem.The experimental and theoretical results showed rea-sonable agreement, indicating that photoelastic tomography can be applied to measuring stress fields in soft materials.Reconstruction of all stress components allows, for example, analysis of the von Mises stress, which is used for the evaluation of the yield value of a material.
To validate the approach for dynamic Hertzian contact, in which a rigid sphere impacts a gel, the dynamic stress field in the gel was reconstructed from the photoelastic parameters measured.Consequently, a wide range of stress fields from O(10 −1 ) to O(10 1 ) kPa were measured simultaneously with large deformation.The current measurement system can also stress wave propagation, in which the stress tensor components change rapidly with time.In addition, the normal impact force acting on the surface was calculated from the reconstructed stress field, and the scaling law for its maximum value was compared with the theoretical prediction for dynamic Hertzian contact.The experimental results showed good agreement with the theoretically predicted scaling laws, confirming the accuracy of high-speed photoelastic tomography.
We believe that this study has provided an optical methodology for measuring dynamic axisymmetric stress fields in soft materials.This could lead to a better understanding of phenomena related to various engineering processes, such as the stress distributions in materials caused not only by the impact of rigid spheres but also by the impact of droplets or liquid jets.

Figure 1 :
Figure 1: The relationship between the three-dimensional stress field in a material and the photoelastic parameters measured by a polarization camera.

Figure 2 :
Figure 2: (a) A general three-dimensional stress field.(b) The considered element ABX, enclosed by the dashed line.

Figure 3 :
Figure 3: A schematic of the algorithm of the onion-peeling method.The relationship between the axisymmetric field p(r) to be reconstructed and the projection data P (r).
rz acting on the N th ring can be determined.The shear stress of the N -1th ring σ (N −1) rz , using the previously determined shear stress of the N th ring σ (N )

Figure 4 :
Figure 4: A schematic of the experimental setup for (a) static Hertzian contact and (b) dynamic Hertzian contact.

Figure 5 :
Figure 5:  The relationship between the surface deformation δzmax of the gel and the loading normal force Fn when the sphere, which has a radius R of 7.35 mm, is pressed against it.The dashed line indicates the theoretical normal force from Eq. (37) using the measured surface deformation δzmax and the Young's modulus of E = 47.4 kPa.

Figure 6 :
Figure 6: An example of the measured (a) retardation and (b) orientation fields.The loading normal force Fn of 0.55 N is induced by the pressing of a sphere with radius R = 7.35 mm.The black parts in the upper middle area indicate the pressed sphere.(c) and (d) are the theoretical retardation and orientation fields, which can be obtained by integrating the theoretical stress field calculated by solving the Hertzian contact problem with the stress-optic coefficient C = 1.14 × 10 −9 Pa −1 .The white parts in the upper area indicate the deformed surface of the gel due to sphere pressing.

Figure 7 :
Figure 7: The stress-optic coefficient C with different loading normal forces Fn.The dashed line indicates the mean value of the stress-optic coefficient of the urethane gel used in this study, which is 1.14 × 10 −9 Pa −1 .

Figure 8 :
Figure 8: A comparison between the (left-hand side) theoretical and (right-hand side) experimental stress fields.(a), (b), (c), and (d) show the shear, axial, radial, and hoop stresses, respectively.(e) The von Mises stress calculated from the stress components via Eq.(40).The loading normal force Fn is 0.55 N, induced by pressing the of a sphere with radius R = 7.35 mm.The white parts in the upper area indicate the deformed surface of the gel due to sphere pressing.

Figure 9 :
Figure 9: The line profiles of the reconstructed stress field compared with that of the expected stress field at z/R = 0.7.(a) The shear and axial stresses.(b) The radial and hoop stresses.(c) The von Mises stress.The loading normal force Fn is 0.55 N, and is induced by the pressing of a sphere with radius R = 7.35 mm.

Figure 10 :
Figure10: The temporal evolution of the retardation and orientation fields during the impact of a plastic sphere.The impact velocity of a sphere with radius R = 2.98 mm is V = 2.2 m/s.The white parts in the middle area indicate the deformed surface of the gel due to the sphere impact.

Figure 11 :
Figure 11: The temporal evolution of the (upper row) shear and (lower row) axial stresses during the impact of a plastic sphere.The impact velocity of the sphere with radius R = 2.98 is V = 2.2 m/s.The white parts in the middle area indicate the deformed surface of the gel due to the sphere impact.

Figure 12 :
Figure 12: The temporal evolution of the (upper row) radial and (lower row) hoop stresses during the impact of a plastic The impact velocity of the sphere with radius R = 2.98 mm is V = 2.2 m/s.The white parts in the middle area indicate the deformed surface of the gel due to the sphere impact.

Figure 13 :
Figure 13: The kymographs of (a,b) the shear and(c,d) the axial stresses acting on the substrate surface, σrz(r, z = 0, t) and σzz(r, z = 0, t), for a sphere radius R = 2.98 mm and impact velocity V = 2.2 m/s, with different color bar limits.

Figure 15 (
Figure15(a)  shows the relationship between F n,max /ρV 2 R 2 and ρV 2 /E with Eq. (49) shown as a dashed line.The trend of the experimental peak force agrees with the theoretical prediction.Additionally, we can obtain the scaling law for the peak time t Fn,max from the relationship z max ∝ t Fn,max V and Eq.(47), as

Figure 14 :
Figure 14: The temporal evolution of the normal force acting on the surface during the impact of a sphere at different impact velocities.

Figure 15 :
Figure 15: (a)  The relationship between the peak force Fn,max nondimensionalized by the inertial force ρV 2 R 2 and ρV 2 /E.The dashed line shows the scaling law for the peak force (Eq.(49)).(b) The relationship between the peak time tF n,max nondimensionalized by the impact time T = R/V and ρV 2 /E.The dashed line shows the scaling law for the peak times (Eq.(51)).