Development of a crushable foam model for human trabecular bone

Finite element (FE) simulations can be used to evaluate the mechanical behavior of human bone and allow for quantitative prediction of press-fit implant fixation. An adequate material model that captures post-yield behavior is essential for a realistic simulation. The crushable foam (CF) model is a constitutive model that has recently been proposed in this regard. Compression tests under uniaxial and confined loading conditions were performed on 59 human trabecular bone specimens. Three essential material parameters were obtained as a function of bone mineral density (BMD) to develop the isotropic CF model. The related constitutive rule was implemented in FE models and the results were compared to the experimental data. The CF model provided an accurate simulation of uniaxial compression tests and the post-yield behavior of the stress-strain was wellmatched with the experimental results. The model was able to reproduce the confined response of the bone up to 15% of strain. This model allows for simulation of the mechanical behavior of the cellular structure of human bone and adequately predicts the post-yield response of trabecular bone, particularly under uniaxial loading conditions. The model can be further improved to simulate bone collapse due to local overload around orthopaedic implants.


Introduction
Trabecular bone is a spongy, porous material with a cellular structure. It is present at the end of all long bones, such as the femur, tibia, and humerus, in the vertebral bodies of the spine, and in flat and irregular bones such as the pelvic bones [1]. Bone mineral density (BMD), age, sex, geometry, and anatomical site all have an influence on the material properties of this type of bone [2]. Mechanical properties of trabecular bone play an essential role in the biomechanical response of the whole skeleton.
The prediction of the biomechanical response of trabecular bone can be used to study the fracture risk or bone collapse, but also the press-fit fixation of orthopaedic implants. The mechanical response of the periprosthetic bone has a direct effect on implant fixation [3]. Already during implantation, plastic deformation of trabecular bone may occur, but also due to local overload of implants, such as the collapse of a tibial tray or unilateral knee prosthesis [4,5]. For the development and design of new devices, it is imperative to understand these nonlinear interactions between bone and implant. Nonlinear finite element (FE) analysis can be used to study the situations mentioned earlier. However, an adequate material model that captures the post-yield behavior of the trabecular bone is necessary for accurate results [6,7].
Several material models have been used to replicate the biomechanical response of human bone. The post-yield behavior in these models mainly has been based either on a softening von Mises (sVM) criterion, a Drucker-Prager (DP) criterion, a Mohr-Coulomb (MC) criterion, or a crushable foam (CF) model [8].
An ideal plastic model based on a sVM yield criterion was developed by Keyak et al. [9] for femoral fracture prediction in case of metastatic lesions. This model was further improved by Keaveny et al. [10] by introducing tension-compression strength asymmetry in the elastic-perfectly-plastic material model. In 2006, a DP model was used to simulate the compression and tension behavior of cortical bone, which was later modified for trabecular bone compression [11]. Also, an extended DP model was calibrated in 2009 by Mullins et al. [12] to capture the post-yield behavior of the bone, considering a pressure-dependent yield effect. An MC criterion was used by Tai et al. [13] and Wang et al. [14] to investigate bone strength based on bovine cortical bone and canine cancellous bone, respectively.
A CF model with isotropic hardening was developed by Deshpande and Fleck [15]; they introduced and calibrated a hardening equation based on hydrostatic and uniaxial compression data, defining an elliptical yield surface for metallic foams. Kelly and McGarry [16] demonstrated the CF model's application for bovine trabecular bone and rigid polyurethane foams [16]. They showed that this type of continuum constitutive model could be implemented for the simulation of bone behavior. The CF model is better able to incorporate the effect of pressure dependency on the yield surface than the DP or MC criterion [16]. The identification of a CF parameters with constant values for human femoral bones and vertebral bodies was performed in a study by Kinzl et al. [17] based on the yield data of whole bones. They captured the mechanical response of the bones until the ultimate yield point. In 2019, Schulze et al. performed a study in which a CF model was calibrated for synthetic bones [18]. It was shown that the CF plasticity model provided accurate predictions of acetabular cup deformation.
To develop the CF models for trabecular bone, one needs to obtain the mechanical response in two loading configurations: uniaxial compression and confined compression. As Zhao et al. reported in 2018 [19], mechanical testing of human bone under uniaxial conditions has been widely investigated over the years by numerous authors. Keyak et al. and Keaveny et al. [6,[20][21][22] performed comprehensive studies on the prediction of post-failure behavior of bone in uniaxial configurations. In the study by Kelly and McGarry [16], the effect of hydrostatic pressure on post-failure behavior of bovine trabecular bone was investigated, which was demonstrated to be of importance under physiological loading conditions. However, there is a limited number of studies that focus on the effect of hydrostatic pressure on post-yield behavior of human trabecular bone. Therefore, the goal of this study was to determine the post-yield response of human trabecular bone under uniaxial and hydrostatic conditions. The present study comprises an experimental and numerical evaluation of the mechanical response of human trabecular bone. In the first phase, through experimental investigations, we determined the required parameters for the CF model (Young's modulus, yield stress, and Poisson's ratio). We examined trabecular bone samples taken from human cadaveric tibias in two configurations, uniaxial and confined compression. In the second phase of the project, the CF model was implemented in FE analyses. A numerical simulation was performed to predict the mechanical behavior (elastic and plastic) of bone under both uniaxial and confined compression. Then, the FE results were assessed in comparison to the experimental data.

Material and methods
This study is divided into three main parts: 1-experimental testing of bone specimens under uniaxial and confined compression, 2-identification of material parameters for the CF model based on the mechanical testing data, and 3-FE simulations based on quantitative computed tomography (QCT) and comparison of the results with the experimental data.

Specimen preparation
Sixty-two cylindrical bone samples taken from the proximal site of 8 fresh-frozen human cadaveric tibias were examined. There were three female and five male bones, and the age range was from 62 to 93 years, with an average of 69 years. The proximal part of all tibias, including the menisci and condyles, were removed parallel to the plateau (Fig. 1a). Cylindrical samples were taken from the trabecular bone sites (Fig. 1b). QCT scans of the drilled bones were captured with a voxel size of 0.4 × 0.2 × 0.2 mm (Toshiba Medical Systems, Tokyo, Japan − 120 kV 260 mA) while a solid calibration phantom (Image Analysis, Columbia, KY) was placed under the samples in the scanner [23].
All samples were milled to the same height of 12 mm and a defined diameter of 11.6 mm, providing an aspect ratio (i.e., height to width ratio) of 1.03, avoiding buckling of the samples during the compression test [19].
To ensure fluids would not interfere with the digital image correlation (DIC) measurements, the bone marrow and other fluids were removed by centrifuging the specimens.

Mechanical setup
The specimens of each bone were divided into two equal groups, which were subjected to either uniaxial or confined compression. Uniaxial compression was performed by placing the specimens between two fixed custom-made parallel platens (Fig. 2b).
The confined compression test was performed using a custom-made chamber and plunger. The inner diameter of the chamber was 11.60 mm and the plunger had an outer diameter of 11.58 mm. The measured diameters of the confined samples varied from 11.25 to 11.6 mm, which resulted in a mean tolerance of 0.10 mm (SD = 0.11). The compressive load was applied through a plunger with a ball joint connection to the load cell (Fig. 2c).

Loading
All specimens were subjected to a preload of 5-10 N to ensure uniform contact between bone and platens. Subsequently, the specimens were subjected to a destructive compression load at a constant velocity of 5 mm/min (strain rate of 0.007 s − 1 ) up to 58% of strain.

Data acquisition
The DIC technique was used to measure bone strains in axial and transverse directions. The crosshead displacement of the testing machine was calibrated against the DIC displacement to account for machine compliance in the setup [24]. The deformations were obtained during the uniaxial compression test based on the DIC data to calculate the Poisson's ratio. Images of the uniaxial compression test were continuously captured and deformations of the samples were calculated based on a custom-written Matlab script (Matlab 7.12.0 (R2018a), Mathworks, MA, USA) using automatic edge detection. It was assumed that the cylindrical bone samples deformed in an axisymmetric manner during the experiment.
The force-displacement data of the compression test was converted to nominal stress-strain curves based on the original dimensions of the specimens. Three significant mechanical properties were measured ( Fig. 3): 1. uniaxial compressive stiffness of the elastic region (Young's modulus, E) 2. the yield point, based on a 0.2% strain offset (σ) (for both uniaxial and confined configurations) and, 3. the ultimate compressive stress. After omitting the initial toe region and the yielding region of the curve [6], the Young's modulus (E) was computed based on the elastic portion of the stress-strain curve. The yield stress (σ) was calculated using a 0.2% offset from the elastic line, and the ultimate stress was defined as the peak of the stress-strain curve within the strain-range lower than 10%. Of these parameters, the Young's modulus, yield stress, and Poisson's ratio were applied to the FE model, while the ultimate stress was measured only for comparison and validation of the computational results.

Material parameter identification
The Crushable Foam model with an isotropic hardening rule (ICF) is governed by the von Mises equivalent stress (q) and the hydrostatic pressure (p). The yield surface is an ellipse centered at the origin in the p-q stress plane. The yield surface extends along the pressure axis under the hydrostatic state (Fig. 4).
The yield surface of the CF model with isotropic hardening (F ICF ) is given by: where B is the size of the q-axis of the yield ellipse, σ uc is the absolute compressive strength under uniaxial loading and a is the shape factor of the yield ellipse and is defined as: In these equations K is the compression yield stress ratio, and p 0 c and σ 0 uc are the initial yield stress under hydrostatic and uniaxial compression, respectively. The shape factor a defines the relative magnitude of the axes of the yield ellipse in the p-q stress plane. In the particular case when a = 0, the crushable foam model corresponds to the von Mises yield criterion. Moreover, the flow potential (G) was chosen as [16,25]: where β characterizes the principal axis lengths of the flow potential ellipse in the p-q stress plane and is correlated by the plastic Poisson's ratio as: These relations define the geometry of the isotropic CF yield criterion in the q-p plane. Additionally, to define the corresponding work hardening slope (H) of the evolving yield stress, the following linear relation was used [15]:  where σ e is the von Mises effective stress and σ is the equivalent stress, and h σ and h p are the slope of the stress versus the logarithmic plastic strain curve under uniaxial and hydrostatic compression (see Appendix for more details).
Using the dimensionlessK parameter and compressive strength as a function of ash density or BMD, the model can be applied to the cellular structure of bone.
The initial yield stress under uniaxial conditions was measured by performing a compression test on the specimens. To obtain the initial hydrostatic yield point, triaxial compression test data is required. Considering Hooke's Law and Poisson's ratio of trabecular bone, the confined principal stresses can be converted to a hydrostatic state, and the yield data can be used to define the K parameter.
Following the method implemented by Keyak et al. [6], the Young's modulus (in uniaxial compression) and yield stress(both for uniaxial and confined compression) of all samples were computed through linear regression of the experimental data. Subsequently, another regression analysis was performed to correlate the calculated model parameters with BMD using Power-law distributions. The obtained correlations were evaluated using Pearson's correlation coefficient (r), and the standard error of the estimate (SEE) was calculated to assess the accuracy of each prediction.

Numerical simulation
In order to assign BMD values to the specimens, a computer model of a cylinder with a diameter of 11.5 mm and a height of 12 mm was created. The cylinder was then virtually placed in the QCT images at the position of the drilled specimens to assign the local Hounsfield units (HU). The obtained values were then converted to BMD using a calibration phantom. The diameter of the cylinder was defined slightly smaller (0.1 mm) than the actual size to reduce the error of positioning.
FE analysis was performed to validate the outcomes of the numerical simulations against the experimental data for five specimens of each configuration. To accurately replicate the heterogeneous BMD distribution of the specimens, 3D models were created for these five specimens. The stress-strain response and yielding patterns of the samples in the simulation were selected to compare with the experimental results. The stress-strain response included compressive stiffness, yield stress, and post-yield behavior. The analysis was done in a displacementcontrolled situation to mimic the experimental condition. A mesh convergence study was performed for a single BMD specimen with four different mesh sizes using 4-noded tetrahedral elements: 0.5, 1.0, 1.5, and 2 mm. The total strain energy of the model was selected as a target measure to determine the differences between mesh sizes. Convergence was assumed for a difference of less than 10% [26], which was achieved for the first three mesh sizes. Considering the mesh convergence results and the optimum mesh size for assigning the BMD values, the cylinder model with an element size of 1 mm was selected. Five specimens from each group with a different BMD, ranging from minimum to maximum values, were selected for the FE simulations, representing the whole BMD range of the specimens.
Simulations were performed in Marc-Mentat (MSC Software Corporation, Santa Ana, CA, USA) with a modified yield criterion for two loading conditions. The analysis was conducted in an implicit scheme with incremental loading and the large strain option was selected as the nonlinear procedure. The bone was considered as a heterogeneous material, and the stiffness values were applied through a compliance matrix for each element. A FORTRAN routine was combined with available user-subroutines of the Marc-Mentat Libraries [27] to define the isotropic CF criterion that was dependent on the BMD values. The algorithm consisted of the definition of (1) the compliance matrix and Hooke's Law (HOOKLW), (2) the yield surface definition (ZERX), (3) the flow potential rule (NASSOC), and (4) the hardening rule (WKSLP).
Computational implementation included a yield stress update algorithm in combination with a hydrostatic pressure dependency (Fig. 5). In this algorithm, first, the material parameters of the constitutive model were imported (ReadF). Next, the compliance matrix of stress-strain was defined based on Hooke's Law (HOOKLW). Then, the ZERX subroutine was used to apply the yield criterion based on the current total stress. In case the yield surface value was smaller than zero (F ICF < 0), the computation continued as elastic. In case of yielding, the algorithm would call a subroutine to define the flow direction for plasticity (NASSOC). Subsequently, the WKSLP was used to prescribe the corresponding hardening slope. Finally, the YIEL user subroutine was used to update the value of yield stress based on the current work hardening and the volume change. The analysis was continued iteratively until the solution converged. This approach allowed for considering the hardening-softening behavior of the trabecular bone structure, which led to a distinctive update of the yield surface.

Experiment
The force-displacement data were collected for 59 of the 62 specimens; no data were obtained for three samples. Two specimens were excluded due to structural failure after centrifuging, and one was crushed under the preconditioning load. Thirty-one bone specimens were tested under uniaxial compression and twenty-eight under confined compression.
Typical diagrams of the stress-strain response of trabecular bone samples in the two configurations are shown in Fig. 6a and b.
The obtained power-law equations for the ICF constitutive model are reported in Table 1. The correlation between mechanical properties of the bone and BMD values have been characterized by the K parameter and its corresponding plastic Poisson's ratio. Additionally, the experimental mean values of the measured parameters are also given in this table.

Numerical simulation
The calculated BMD values ranged from 26 mg/cc to 207 mg/cc for the uniaxial samples and 27 mg/cc to 195 mg/cc for the confined samples. The experimental data results, coupled with the numerical simulation outcomes of five selected samples are shown in Figs. 8a, b and 9a, b for uniaxial and confined configuration, respectively.
In the uniaxial configuration, the numerically obtained stress-strain curves were very similar to the experimental data. The model accurately replicated the post-yield trend of the specimens in uniaxial compression. In the confined configuration, as can be seen in Fig. 9, the first part of the stress-strain curve, including the yield and ultimate stress, was consistently simulated in the FE analysis. For three specimens, the experimental results showed a plateau region in the stressstrain curve after 15% strain and started to increase, while in the  other two curves, the ultimate stress was followed by a drop and then an increase in stress. The numerical simulations could capture the stress drops of the two specified samples; however, the amount of this drop was underestimated compared to the experimental results.
Numerically determined, the ultimate compressive stress of the uniaxial and confined configuration was highly accurate compared to the experimental values ( Fig. 10a and b). Note that the ultimate stress of the experimental data was not included as input for the CF material model and was computed based on the FE algorithm as a post-yield response.
The distribution of the equivalent plastic strain (EPS) in the uniaxial configuration is illustrated in Fig. 11 as a representation of plastic yielding. This computed distribution of plastic strain in the FE analysis compares relatively well with the deformations observed in the experiments.

Discussion
In this study, we assessed the mechanical response of human trabecular bone under compressive loading conditions through experimental and numerical simulations. The load was applied in uniaxial and confined configurations to characterize the parameters of a constitutive crushable foam model that was dependent on BMD. QCT-based FE simulations of trabecular bone specimens were compared against experimental results.
The Young's modulus and yield stress measured in the uniaxial configuration of the tibial trabecular bone varied from 31 MPa to 458 MPa and 0.53 MPa to 7.59 MPa, respectively. These values are comparable to Young's modulus (8-1310 MPa) and yield stress (0.83-24 MPa) reported in previous studies [20,[28][29][30]. The measured values of the yield stress under confined compression in the current study (0. 16-9.02 MPa) were within the range of 0.24-31.59 MPa reported by Carter et al. [31]and Charleroi et al. [32].   [15].
Regression analysis showed a strong correlation (p < 0.001) of three measured parameters (compressive stiffness, uniaxial yield stress, and confined yield stress) with BMD values. The values derived here (Table 1) had a good agreement with reported values of the mechanical properties of proximal tibial specimens in the review study by Keaveny et al. [1]. However, our results were quite different from those of Keyak et al. [6], who found a Young modulus ranging from 135 MPa to 1200 MPa, and yield stress ranging from 1.36 to 9.8 MPa. This difference may   be due to several factors. First, the average stiffness of the cadaveric bone samples in the study by Keyak et al. [6] was much higher compared to our samples, which may lead to different regression results. The second factor may be related to differences in the experimental setup. In the study of Ketak et al., the displacement was applied at a rate of 9 mm/min (strain rate of 0.01 s − 1 ) on wet bone samples (including bone marrow), while in the current study, the displacement velocity was 5 mm/min (strain rate of 0.007 s − 1 ), and liquids (bone marrow and fat) were removed from the samples before testing. However, it has been stated that a strain rate of less than 10 s − 1 should not affect the measured stiffness of the samples, with or without the bone marrow [31,33].
Applying the material models to FE simulations of trabecular bone in an accurate manner is a challenging procedure [8]. Many attempts have been made to apply constitutive continuum models to trabecular bone, including softening Von Mises (sVM), Drucker-Prager (DP), Extended Drucker-Prager (EDP), Mohr-Coulomb (MC), and Crushable Foam (CF) models [9][10][11][12][13][14]16,17]. For constitutive modeling of cellular structure, it has been stated that the hydrostatic and deviatoric stresses should be considered together in FE analysis to have an accurate response [8,16]. However, among the models mentioned above, only the CF and EDP are functions of hydrostatic and deviatoric stress.
The present study is the first attempt to characterize the material parameters of the ICF model dependent on the BMD. Several attempts have been made previously to calibrate the CF model for human [7,17] and for synthetic [18,34] bone, using a single value for the K parameter. A simplified form of a calibrated CF model (with a constant value for the corresponding parameters) was applied to femoral bones and vertebral bodies by Kinzl et al. [17]. They showed accurate results for the prediction of ultimate strength, as well as the damage distribution. In contrast to the current study, they did not include hardening-softening equations in their CF model. Therefore the computed load-displacement curves in their study were only valid until the ultimate force. The CF model that was used in the study of Kelly et al. [7] to investigate the behavior of vertebral bodies showed a good agreement with experimental results. However, the μCT model in their study did not allow for using a realistic 3D geometry in a macroscale analysis. The values for the K parameter reported in these previous studies were in the range of 0.85 to 1.33, which is comparable with the values in the current study.
As shown in Fig. 8, using the ICF model, the mechanical behavior of trabecular bone was adequately reproduced in numerical simulations and captured the crucial points of the stress-strain curve under compression. The measured initial yield stress and Young's modulus were already applied to the FE model as input data and had a strong correlation with the numerical outcomes. Interestingly, although no parameters from the post-yield region of the experiment were considered as input data, the numerical simulation could accurately replicate the plastic behavior of the specimens (Fig. 8b). The post-yield region of the stress-strain response (ultimate stress point, softening part, plateau region, and the hardening part) could be quantified well in the simulation results. Given the three mechanical parameters as the initial input for the simulations, the ultimate stress was computed in FE simulations, which were very similar to the experimental results (Fig. 10).
Considering the confined compression, Kelly and McGarry [16] showed a decreasing stress after the peak point in their bovine CF model, while their experimental results showed an increasing stress after the peak point. In the current study, the numerical outcomes of the confined simulations showed a good agreement with the experimental results up to 15% strain (Fig. 9). After this strain level, the computed stress increased relative to experimental values due to the overestimated hardening rule. According to the study of Yu et al. [35], for a realistic FE simulation of confined compression, it is necessary to obtain sufficient modifications for the constitutive formulation. They stated the hardening-softening rule must be dependent on the confined pressure, and the flow rule must be dependent on confinement level and incremental rate. Therefore, in order to apply these modifications to a FE analysis, a considerable amount of experimental data under different levels of confinement is required. If these data are not included in the plasticity model, the simulations of confined configurations lead to an overestimation of the hardening rule (as is the case in the current study) or softening rule (e.g. the study of Kelly and McGrey [16]).
Since the material properties in the models were based on the actual BMD, the localized plastic behavior resulted in a good prediction of the yielding pattern. The deformations of the experimental specimens under uniaxial compression were qualitatively similar to the plastic strain distributions seen in the simulations (Fig. 11). By considering the flow potential rule of the CF plasticity model (Eq. (5)), the direction of the plastic strain rate vector was updated independently of the yield surface. Describing the differential changes of the plastic strain component based on the flow potential rule allowed for realistic deformations of the trabecular bone. Therefore, identifying the plastic zones with the maximum values of the equivalent plastic strain in numerical simulations could indicate the failure pattern of the experimental specimens.
Kelly and McGarry [16] applied the DP, MC, sVM, and CF model to simulate the compression situation of bovine trabecular bone. They demonstrated that none of these constitutive models could capture the confined compression response of trabecular bone, except for the CF model. The parameters K and v p in their study were calibrated against the average stress-strain curve of the experimental results. Therefore, these parameters are only dependent on one specific BMD value of bovine bone and not on the entire range. In the study of Schulze et al. [18], it was shown that the CF plasticity model with an isotropic hardening rule resulted in an accurate prediction of deformations in synthetic bone. Considering the fact that the CF model mainly depends on the K parameter (the ratio of uniaxial yield over hydrostatic yield), they calibrated their model by assuming equality of all the yield ratios and set this parameter equal to 1. However, the current study shows that a variable value of the K parameter results in a more accurate response of the cellular structure of trabecular bone.
Although this study shows good agreement between the experiments and the numerical simulations, there are some limitations. Minor experimental errors are difficult to avoid. These errors may be related to the end effects of the platens, machine compliance, sample preparation, structural damage of specimen, and the aspect ratio of the specimen's geometry, which have been reported previously as standard experimental errors [19]. Also, the perfect alignment of the bone samples with the cutting direction and axial loading remains challenging, due to the complex structure of the trabecular bone. Regarding the numerical simulation, in the uniaxial configuration, the model could accurately predict the mechanical response of the trabecular bone. However, in the confined configuration, the stress was overestimated after 15% strain. The isotropic hardening rule of the CF model contained pressure-dependent parameters. However, according to Yu et al. [35], to simulate pure confinement in FE analysis, the hardening-softening rule should be dependent on the confining pressure and requires sufficient experimental input, which was lacking in the current study. Further experimental tests and simulations of the confined situation are required to further improve the CF model. Although it is fundamental to perform hydrostatic tests (or confined tests in this study) to obtain the CF parameters, the confinement boundary condition in this level is not necessarily a correct representation of the physiological conditions in human bone. In reality, collapsing bone is surrounded by other deformable bone (both trabecular and cortical bone). Hence, a realistic mechanical response of the post-yield behavior (and its modeling) will probably be somewhere between the uniaxial and confined conditions where the confining level depends on the stiffness properties of the surrounding bone.

Conclusion
In the present study, mechanical properties of human trabecular bone were experimentally determined dependent on BMD. Using these properties, an isotropic crushable foam model was developed. This model realistically predicted the post-yield behavior of trabecular bone under uniaxial conditions. Also, it could reproduce the pure confined compression until 15% of strain. The CF model can properly simulate orthopaedic device performance, particularly focusing on bone collapse due to the local overload around orthopaedic implants.

Ethical approval
Not required.

Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of Competing Interest
Authors declare that they have no conflict of interest.
where c and d are experimental constants, ρ BMD is the density of an element, M is the total mass and V 0 is the primary volume. The relation between ε p and ε p vol can be derived based on the reported equations in Deshpande and Fleck [15], by defining the ratio of q over p as R: (9) σ uc = σ 0 consequently, the following is obtained: ∂σ uc ∂ε p = − α 2 cd(ρ BMD0 ) d (R + α 2 ε p ) d+1 (11) in which ε p is calculated each increment, and other parameters are given as input data. The same procedure was carried out for h p by converting the confined compression stress (σ cc ) to the hydrostatic stress (p) using Hooke's Law for isotropic material: