40

Using the deﬁnition of the Poisson’s ratio (ν), we can relate the transverse strains with the axial strain ε11 using the relation presented in Eq. 3.1. (1 + ε22 ) = (1 + ε33 ) = (1 + ε11 )−ν ,

(3.1)

where the principal stretches are deﬁned in terms of the principal nominal strains by: λi = 1 + εii . Additionally, the deformation gradient tensor (F ) can be also expressed in terms of λi as shown by Eq. 3.2. To satisfy the assumption of an incompressible material, the Jacobian (J) should be equal to 1. Therefore, J = det(F ) = λ1 λ2 λ3 = 1, where λ1 0 0 F = 0 λ2 0 . (3.2) 0 0 λ3 Considering that this procedure will be the same for the rest of the material models, let us ﬁnd the analytical solution for a Neo-Hookean material which is the simplest to be determined. From Eq.2.31 and Eq.2.39 we can deduce the stress-strain relationship from the strain energy density of the Neo-Hookean material model [6]: σij =

C10 1 (Bij − Bkk δij ) + K1 (J − 1)δij . 5/3 3 J

(3.3)

Now using Eq. 2.6 and Eq. 3.2, the normal stress for an incompressible Neo-Hookean material is given by Eq. 3.4. 1 σ11 = 2C10 (λ21 − (λ21 + λ22 + λ23 )), (3.4) 3 −1/2

(based on Eq.3.1). where λ2 = λ3 = λ1 Using the relationship between strains and stresses during a uniaxial compression test, the normal components of the Cauchy Stress can be deﬁned in terms of the principal stretch in direction of the load application, as shown in Eq. 3.5 and Eq. 3.6. 4 C10 (λ21 − λ−1 1 ), 3

σ11 =

(3.5)

2 2 C10 (λ−1 (3.6) 1 − λ1 ). 3 Now, knowing the relation of the Cauchy Stress with the Nominal Stress (see Eq. 2.16) we can deﬁne the following nominal stresses: σ22 = σ33 =

S11 =

4 C10 (λ1 − λ−2 1 ), 3

(3.7)

2 C10 (λ−2 (3.8) 1 − λ1 ). 3 Finally, and using Eq. 2.6, Eq. 3.7 and Eq. 3.8, the constitutive equation (analytical solution) of an incompressible Neo-Hookean Material submitted to an uniaxial compression test is given in terms of the Nominal Stress diﬀerences by Eq. 3.9. We rather use Nominal Stress instead of the Cauchy Stress, because Sij corresponds to the internal force per unit of undeformed area acting within a solid, which is easier to determine during an experimental procedure. The constitutive equations for the Reduced Polynomial, Mooney Rivlin and Ogden models are similarly found and determined by Eq. 3.10, Eq. 3.11 and Eq. 3.12, respectively. S22 = S33 =

Neo-Hookean: nd

2

S11 = 2C10 (λ1 − λ−2 1 ),

Order Reduced Pol.:

S11 =

Mooney-Rivlin Solid:

S11 =

Ogden Form:

S11

¯ 2(λ1 − λ−2 1 )(C10 + 2C20 (I1 2(1 − λ−3 1 )(C10 λ1 + C01 ),

µ1 1 −1 1 −1 = 2 (λα − λ−0.5α ). 1 α1 1

(3.9) − 3)),

(3.10) (3.11) (3.12)

CHAPTER 3. CHARACTERIZATION OF SOFT TISSUE

41

Experimental Setup: Uniaxial Compression Test As shown by Francheschini et al. [24], human brain tissue deforms similar to ﬁlled elastomers, and it should be modeled as a nonlinear solid with small volumetric compressibility. Girnary [29] also highlight that silicone brain phantoms provide good results to simulate brain behavior. In the present study, we determined the parameters for four diﬀerent hyperelastic models to simulate the mechanical behavior of a platinum-cure silicone rubber submitted to a compression test. The silicone rubber is called Ecoﬂex 00-10 (from Smooth-On, Inc.) and with the aim of having a closer behavior to brain tissue, we used a proportion (20%) of a softener called ”Slacker” (also from Smooth-On, Inc.). We used a tissue phantom rather than biological tissue because this facilitates obtaining repeatable results in a controlled environment. Each component of the rubber solution was evenly mixed according to the recommendations of the manufacturer and then formed in a cylindrical mold of 80 mm diameter and 70 mm height. We selected this dimensions based on the size of the “Truth Cube”. A Force/Torque sensing system was installed in a laparoscopic grasping tool and was calibrated under a rigorous accuracy testing by the manufacturer (ATI, Industrial Automation). The certiﬁcate of calibration establishes that the maximum amount of error for the Z axis (the needle axis) is 0.75% of its full-scale load. This calibration was done at 74◦ F (22◦ C) and the experiments of the entire thesis were done at this temperature. We used a force/torque sensor, (ATI Mini40 SI-40-2) with 0.02 N resolution, attached to a laparoscopic grasping tool, which was ﬁxed to a rigid plate (see Fig. 3.4). The contact areas were lubricated with talcum powder in order to minimize lateral friction. The plate was displaced using a Phidget Bipolar Stepper Motor to compress the tissue at a constant velocity of 0.248 mm/s, until the plate was displaced by 10 mm. We selected an stepper motor which provided a controlled displacement for low velocities and its maximum speed was 0.248 mm/s. See Appendix C for more information in the parameters of the motor and F/T sensor. Based on the studies of brain tissue from the University of Western Australia (leaded by Professor Karol Miller [52], [54], [51]), the research done by Franceschini [24] from the University of Trento, and Kohandel [41] from the University of Waterloo, we determined the velocity of our experiments which give us low strain rates typical for neurosurgical procedures. Methodology We created a program in C++ that integrates the measurements of Forces (given by the F/T sensor) with Displacements (given by the motor) using the time (in milliseconds) as a reference. We assumed the experiment occurs under symmetric conditions, with a homogeneous, isotropic and incompressible material, and that imposed deformations were small compared with the original size of the cylinder. We selected the hyperelastic models deﬁned by the equations 3.9, 3.10, 3.11 and 3.12 to estimate the mechanical response of our specimen. To calibrate the material parameters, we minimized the square of the absolute error between the stress-strain curve obtained with the analytical solution with respect to the experimental measurements (see Eq. 3.13 and Fig. 3.5). This estimation is performed using a least-squares ﬁt. We used absolute errors instead of relative errors, because they provide a better ﬁt for large strains (deformations bigger than 0.05). m ∑ min{absolute error2 } = min{ (Sexperiment − Sanalytical )2 },

(3.13)

i

where m is the total number of data points. The experimental force displacement relationship was found with the integration of the measurements given by the motor and the force sensor. The experimental Nominal Stress (S11 ) was found by dividing the reaction force at every sample point by the undeformed contact area. The

CHAPTER 3. CHARACTERIZATION OF SOFT TISSUE

42

Figure 3.4: Experimental setup for the compression test of a cylindrical shape of silicon rubber (Ecoﬂex -0010). The surfaces were lubricated and the tissue was compressed to a strain of 0.14 Nominal Strain corresponded to the displacement of the plate divided by the undeformed height of the cylinder. Once we obtained the material parameters, we evaluated its stability using the Drucker Stability Test, because the material should satisfy the fundamental assumptions underlying continuum constitutive equations [6]. This stability criterion for an incompressible material establishes that the work done by the tractions through the displacements is positive or zero. This condition is satisﬁed when the stress-strain relation is bigger or equal to zero (see Eq. 3.14). ∆τij ∆εij ≥ 0,

(3.14)

where ∆τij is the change in Kirchhoﬀ stress, and ∆εij is an inﬁnitesimal change in displacement. “The stability criterion is violated wherever the stress decreases with strain in tension, or increases with strain in compression” [6][66]. For the Neo-Hookean material model, the Drucker stability is satisﬁed when the coeﬃcient C10 is positive. Results The stress-strain experimental curve for the uniaxial compression test and the four predicted curves using the diﬀerent hyperelastic models are shown in Fig. 3.6. To calculate the error, we picked data points at evenly spaced strain intervals over the range of strains that we will require later to

CHAPTER 3. CHARACTERIZATION OF SOFT TISSUE

43

Figure 3.5: Method to characterize diﬀerent material models based on experimental data and the analytical solution od the problem.

Figure 3.6: Stress-Strain relationship for a standard compression test and its comparison with the analytical solution of the problem using diﬀerent hyperelastic material models. do needle indentations. This gave us similar ﬁtting accuracy throughout the entire strain range. To facilitate the interpretation of the results, we also reported in Table 3.3 the R-Squared values for each ﬁtting result with their corresponding parameters (we organized them from best to worst). The Drucker Stability test showed that the four models were stable for all strains. Table 3.4 shows the corresponding initial Young’s Modulus (E0 ) and initial Shear Modulus (µ0 ) obtained with each of the ﬁtting parameters (using the relations given by equations 2.36, 2.38, 2.39

CHAPTER 3. CHARACTERIZATION OF SOFT TISSUE

44

Material Model

Parameters [Pa]

R-Squared

Reduced Polynomial (N=2)

C10 = 1404.9, C20 = 699.4

0.99971

Mooney-Rivlin

C10 = 624.3, C01 = 746.1

0.99897

Ogden

µ1 = 3225.9, α1 = 14.8

0.99888

Neo-Hookean

C10 = 1461.6

0.99875

Table 3.3: Material parameters and the corresponding error for diﬀerent material models using a standard compression test measurements and 2.41). Material Model

E0 [Pa]

µ0

Reduced Polynomial (N=2)

8429.6

2809.8

Mooney-Rivlin

8223.0

2741.0

Ogden

9677.9

3225.9

Neo-Hookean

8769.6

2923.2

Table 3.4: Initial Young’s and Shear modulus obtained with a Standard Compression test and the analytical solution of the problem.

Summary In this section we obtained the material properties of a cylinder of Silicone Rubber which contained 20% of softener (Slacker) using the analytical solution of a standard compression test. We used four material models (Neo-Hookean, second order Reduced Polynomial, Mooney Rivlin and Ogden Form) and we concluded that the Reduced Polynomial gave the best ﬁt (R − squared = 0.99971) at the same time of being stable based on the Drucker criterium. The parameters for this model were: C10 = 1404.9 and C20 = 699.4 which lead us to an initial Young’s Modulus of 8.4 kP a. We conclude that the use of the analytical solution is only feasible for very simple approaches (like the one we just presented in this section), because the solution of the diﬀerential equation that relates stress with strain is not simple to be obtained for all tool-tissue interactions. In the following section, we will compare and validate the obtained material properties with other calibrations using Inverse FEM over the same specimen.

3.2.2

Characterization of soft tissue using Inverse FEM and a bonded compression test

In complex scenarios, the use of an analytical solution in the calibration of a material is not feasible. Therefore, solving the stress-strain diﬀerential equation by means of a numerical approximation emerges as a good alternative. This section includes the calibration of the same material of section 3.2.1 but using bonded compression test experimental data and Inverse FEM. We validated the parameters by comparing them with the results found in the previous section. Inverse Finite Element Method As mentioned in section 2.3, the FEM give us a numerical approximation to the solution of a partial diﬀerential equation. Usually to solve this type of problems, we know the boundary conditions, material properties and geometry before we obtain the stress-strain (or force-displacement) at any node in the domain. Now let us assume we do not know the material properties, but we know the

CHAPTER 3. CHARACTERIZATION OF SOFT TISSUE

45

force-displacement relationship (based on experimental measurements), geometry and boundary conditions. The Inverse Finite Element Method, allow us to ﬁnd the material properties by running multiple FEM simulations, until we get the optimal force-displacement relationship that is closed to the experimental measurements. This procedure is showed in Figure 3.7.

Figure 3.7: Inverse Finite Element Method to characterize diﬀerent material models based on experimental data.

Bonded Compression Test The use of compression testing to determine the mechanical behavior of a material is widely used. In contradiction to a simple compression test, the bonded compression test does not include the use of lubricant between the rigid plates and the phantom tissue. As a result the rubber is compressed and it expands laterally; this is called the “barrelling eﬀect” (see Figure 3.8.) Previous works have determined the analytical solution of this problem for an Elastic material [15], [28], [77]. However, the solution for more complicated material models (such as viscoelastic and hyperelastic) is not trivial.

Experimental Setup: Using the same specimen than the previous section (a cylinder made of Ecoﬂex 00:10 with 20% of softener), we applied compressive forces without the use of lubricant and we executed the experiment three times to validate the measurements. The force sensor, motor, temperature, displacement of the plate, and the velocity of movement, were the same than in section 3.2.1. Method Based on the resultant experimental measurements (speciﬁcally, the force/displacement curves) in conjunction with a FEM simulation, we determined the material properties using the method described in Fig. 3.7. We ran multiple FEM simulations keeping the same geometries and boundary conditions than in the previous simulation, with the diﬀerence than in this case we considered friction between the plates and the tissue. Additionally, we changed the material properties to minimize the error between experimental measurements and the FEM results. From the standard compression test, we deﬁned the material model that better ﬁt the behavior of our silicone rubber, which corresponds to the Reduced Polynomial Model (with N = 2). Therefore, we used this model in the following simulations. Considering that we evenly mixed the

CHAPTER 3. CHARACTERIZATION OF SOFT TISSUE

46

Figure 3.8: Barreling Eﬀect during a Bonded Compression Test. components of the silicone, we assumed the experiment occurs with a homogeneous and isotropic material. Furthermore, we controlled the location of the plates respect to the specimen to ensure symmetric conditions. Because we used a cylindrical shape, we could also analyze this tridimensional problem with an axisymmetric formulation that simpliﬁes the study in a 2D scheme. Based on [20], [21], [65] we determined the friction coeﬃcient between silicone rubber and steel. They deﬁned that this coeﬃcient ranges between 0.6 to 0.9, so we run multiple simulations using friction coeﬃcients of 0.6, 0.7 and 0.9, and we did not ﬁnd signiﬁcant diﬀerences either in the prediction of the reaction forces nor in the geometrical changes. In consecuence, we used a friction coeﬃcient of 0.7 in all our simulations. We simulated bonded compression in ABAQUS 6.10 EF2 as the contact between a rigid plate, the ﬂoor (also an 1D discrete rigid surface) and, an isotropic deformable silicone rubber. The discrete rigid bodies were two lines composed by 10 linear line elements of type RAX2 and 11 nodes. Each rigid body had a reference point (RP) where we applied boundary conditions (BC). The axisymmetric deformable model corresponded to a square of 40 mm by 70 mm, with 500 linear quadrilateral hybrid elements of type C4X4RH and 546 nodes. The mesh size was graded to be more reﬁned close to the areas where we expected higher deformations and coarse near to the axis of revolution (see Fig. 3.9). In all simulations, we allowed nonlinearities, for the material model and for large geometric deformation. We delimited position constraints for the three parts that were provided for our simulation. The edges of the rigid parts were coincident with the upper and lower surface of the deformable body. In the contact we deﬁned a penalty method, where the tangential behavior had the friction coeﬃcient previously determined (0.7). Figure 3.10 displays the boundary conditions applied in this simulation. In blue, one can see the two reference points (RP) for each rigid object and the boundary conditions (BC) are in green. This BC completely restricted the displacement and rotation of the RP-B. The RP-A was displaced 10 mm in y-direction during 40.32 s (which ensured a compression velocity of 0.248 mm/s), its displacement in X and rotation around Z was ﬁxed to be zero. Finally, we also applied a symmetry BC in the X direction. Because this study required a nonlinear analysis, we needed to solve multiple systems of equations. Hence, the solution was found by specifying the loading as a function of time and incrementing time, which allowed us to obtain the nonlinear response. Therefore, for each time increment ABAQUS solved the system of equations using the Newton Method [68]. As for the time integration we deﬁned 250 as the maximum number of increments, 0.00125 was the initial

CHAPTER 3. CHARACTERIZATION OF SOFT TISSUE

47

Y Y

Z

X Z

X

Figure 3.9: Mesh of the axisymmetric model for the FEM simulation of a Bonded Compression Test.

Figure 3.10: Boundary Conditions for the axisymmetric FEM simulation of a bonded compression test increment size, 1.25×10−5 was the minimum increment size, and 0.25 was the maximum increment. Once the solution for a speciﬁc set the material properties was found, an optimization algorithm determined a new set of parameters to redo the FEM simulation. Taking into account the time that each FEM simulation required, we manually deﬁned an initial set of parameters that gave us appropriated results by means of a “Manual Optimization”. With the initial guess (resultant from the manual optimization) we used the Levenberg-Marquardt Optimization algorithm to ﬁnd the material coeﬃcients that better ﬁt the experimental measurements. We used Python to manage inputs and outputs directly from ABAQUS, in order to estimate the error, process the output, modify the material properties, and submit a new simulation until the convergence was achieved. Finally, we validated the stability of the material with the Drucker Stability test (as done in the