Valutazione dei parametri geometrici che influenzano la caratterizzazione poroelastica del tessuto osseo corticale mediante modellazione numerica di prove di nanoindentazione

1. Introduction Nowadays, the mechanical characterization of tissues and other biological materials is of paramount relevance in both clinical and biomaterials science. The mechanical properties have to be analyzed and their links with probing technique length scale and time is topical. In the present work cortical bone is considered as a paradigm of a complex tissue with mechanical functions dependent to the length scale and time. Bone is a biological material composed mainly of collagen fibers and hydroxyapatite crystals. The consituents' arrangement in a precise hierarchical organization contributes to determine the mechanical properties at the macroscopic level. From a structural viewpoint, articular cartilage can be seen as a porous solid saturated with fluid that can flow through pores. This time dependent behavior is due to two main mechanisms: i) viscoelasticity that is intrinsically dependent on the molecular structure and, ii) poroelasticity that is an extrinsic mechanism related to the fluid flow trough pores. The aim of this work is to study the mechanical behavior of articular cartilage with specific regard to the time response and the effects of probing technique length scale. Nanoindentation is an experimental technique that allows investigating the mechanical behavior of a sample in different directions. An indentation test consists in pressing an indenter onto the surface of interest; the tip of the indenter is made of a hard material, usually diamond. The load P and the displacement h of the tip into the material are continuously recorded during a test. The test has two distinct steps, namely the loading and unloading phases, where the maximum load or the maximum displacement are reached before the tip returs to its original position. From the P-h curve it is possible to extract some relevant parameters directly related to the mechanical behavior of the material, as the indentation modulus and hardness. The main advantage of this technique is the absence of stringent requirements in terms of size and shape of specimens, even if the experiments are carried out at nano or micro scales and below 300mN loads. Therefore, instrumented indentation represents a powerful tool for measuring biological tissues properties at the micro and sub-micro levels. While viscoelasticity has been widely studied and closed form solutions for the interpretation of data are available, the interpretation of the poroelastic response still remain an open issue. It’s important to notice the lack of experiments in which the condition of the specimen described in this work are considered and the extremely low availability of FEM simulations where the poro-elastic behavior of cortical bone can be regarded together with its dependence upon the simulation parameters. In addition, it’s necessary to stress how, although there are many trials with specimens in conditions near to the physiological ones, experiments performed with spherical indenters in place of the Berkovich one are not traceable in literature. This latter indenter is not able, as we are going to demonstrate here, to identify the important role of the poro-elasticity in the definition of the mechanical property of bone. In the present work, the time-dependent response of cortical bone, during relaxation and creep experiments carried out with conical (Berkovich) and spherical indenters of several nominal radius, is analyzed. This approach represents an incisive solution because allows the estimation of the local properties and heterogeneities and determination of mechanical response due to the geometry of the indenters, loading conditions and duration of the experiment. The present work aims to simulate experimental nanoindentation tests by adopting a finite element numerical model developed with the commercial code ABAQUS able to capture the main aspects of mechanical poro-elasticity of cortical bone and to adequately reproduce the shape of loading and unloading curve. Thus a comparison between the results of numerical simulations and experimental tests carried out by Carnelli et al. (2011), Wang et al.(2006), Bembey et al. (2006), is presented in terms of agreement between the indentation curves and the outcome values of the investigated parameters, namely indentation modulus and permeability. 2. Materials and methods Poro-elasticity of bone tissue has been studied widely by simulating the nanoindentation tests in different directions. For the material behavior transverse isotropy was assumed for a reference system with the vertical axis coincident with the main axis of the long bone, the radial-circumferential plane is considered as the plane of isotropy. Consistent with the chosen material symmetry, nanoindentation tests performed along the axial direction were simulated using an axisymmetric finite element model, while a three-dimensional model was required to simulate the transverse indentations. Since data from the literature reveal marked differences between yield stress values in tension and in compression for bone, the use of a pressure-dependent yield criterion as the Drucker- Prager one is appropriate. For the numerical simulations, relaxation and creep experiments have been developed for calculating indentation moduli and for identifying the best indenter geometry and the best indenter dimension in order to evaluate the poro-elastic phenomenon and the bone permeability at this scale. These trials have been performed with the computational equivalents of the Berkovich indenter and with 5,7,10,15,21 and 30μm radius spherical indenters. The length of loading step was fixed at 0.1s, the length of the holding one was 10s and the unloading step was firstly fixed at 1s and than at 10s. For relaxation experiments, indentation moduli were extracted by interpolating data from the 90-40% of unloading curve to the Oliver and Pharr (1992) function F= a(h-hf)^m. On the other hand, time-dependent behavior was analyzed by interpolating of creep curves with the analytical solution for a poro-elastic spherical indentation problem defined by Agbezuge and Deresiewicz (1974): H^*=A-A/[1+(T^*/T_0 )^P ] where displacement and time were opportunely normalized: T^*=√ct/(a(t)) H^*=(h(t)-h(0))/(h(∞)-h(0)) where c is a constant parameter depending on permeability and shear modulus. Concerning the elastic properties of the material, elastic constants were extracted directly from the experimental trials of Carnelli et al. (2011) the analytical values of the model were then obtained by introducing the explicit formulation of Delafargue and Ulm (2004). 3. Results The indentation moduli calculated in both the drained (the response is determined by the only solid matrix) and undrained (the response is determined by the fluid component) condition of specimen are reported in Tab 1 where is clear that the ideal and rounded off tip are not able to evidence any remarkable difference between the two condition while the spherical tips do. Indenter M (GPa) Drained Undrained Ideal 27.5 29 Rounded off 25.5 26 R=5μm 23.6 29 R=15μm 23.2 29 R=21μm 23.17 28.7 R=30μm 24 28.5 Tab 1: Indentation moduli extracted from the simulations in both drained and undrained conditions with ideal and rounded off version of Berkovich indenter and spherical indenters used here. The anisotropy ratios extracted from transverse indentation with 15 and 21μm radius tips are consistent with the experimental data reported by Carnelli et. al(2011) as shown in Tab 2: Radius(μm) MT (GPa) MA (GPa) MT/MA HT (GPa) HA (GPa) HT/HA 15 19 23.6 0.8 0.51 0.57 0.89 21 18 23.2 0.78 0.46 0.47 0.97 Tab 2: Indentation moduli in transverse MT and axial MA direction, anisotropy ratio MT/MA, hardness in transverse HT and axial HA direction and anisotropy ratio connected HT/HA for both the indenters here used. The poro-elastic behavior of cortical bone comes out from the creep curves reported. The difference between these curves and the same one derived from a purely visco-elastic material lie in the lack of perfect overlap after time normalization. There was a good agreement between the model and shape of the simulation data for all the creep curve as already noticed by Oyen et al. (2007) for the experimental data of spherical indentation reported by Bembey et al. (2006) so it’s been verified the applicability of the analytical solution proposed by Agbezuge and Deresiewicz for both cases. From the figure shown below it’s possible to observe the development of time constant T0 of the fitting function with the contact radius associated to the different spherical indenters. It stands to reason that, in order to study the poro-elastic behavior of bone and evaluate the permeability at this scale, it’s necessary to choose the indenter’s radius between the point 1 and 2. Then, the permeability value has been studied with creep experiment carried out with a spherical indenter of 7μm. As observable in figure an excessive decrease of the permeability value doesn’t cause any change in data results. These considerations are consistent with experimental and theoretical data reported in section 2.4 and then a permeability value less than O(10-25) isn’t faithful for cortical bone. Even a value ten times bigger than the first one used in these simulations is wrong for the tissue at this scale because the poro-elastic phenomenon associated is too rapid. Thus, the intrinsic permeability value acceptable for cortical bone is included between O(10-23) e O(10-25) m2. 4. Conclusions The results obtained in the present work lead to an important conclusion: if poroelasticity is one of the main dissipation mechanisms for cortical bone exhibiting a time dependent response, then the results in terms of normalized time will be independent of constitutive parameters. Consequently, if further size dependence is still observed in the dimensionless time domain, then further dissipative mechanisms are playing a role and intrinsic viscoelasticity would probably be involved. In this thesis work a method for the interpretation of poroelastic response of cartilage has been introduced, by using nanoindentation tests with spherical tip. This mechanical test is able to show and discriminate viscoelastic and poro-elastic behavior of a complex biological tissue as cortical bone and to extract constitutive parameters by applying numerical and analytical models. In this work we had specific regards for permeability, indentation modulus and anisotropy ratio. The utility of this work lies on: the introduction of a finite element model useful to emphasize the important role of poro-elastic component in mechanical time-response of cortical bone, the identification of the best geometry and dimension of indenter that could understand the poro-elastic phenomenon and the institution of a guide line for the development of new particular experimental protocols able to find and differentiate all the components of bone tissue’s time-dependent behavior.

1. Introduzione Negli ultimi anni sempre più interesse è stato riposto nella caratterizzazione meccanica di tessuti e materiali biologici, sia in ambito clinico sia nello studio dei biomateriali. Per uno studio completo le proprietà del materiale devono essere investigate in correlazione con le lunghezze e i tempi caratteristici delle procedure di prova. Nel presente lavoro l’osso corticale viene presentato come esempio di tessuto complesso le cui proprietà meccaniche sono dipendenti dal tempo e dalle lunghezze geometriche caratteristiche della prova. L'osso è un materiale composto essenzialmente da collagene e cristalli d’idrossiapatite disposti secondo una precisa organizzazione gerarchica che contribuisce a definirne le proprietà meccaniche a livello macroscopico. Da un punto di vista strutturale l’osso corticale articolare è costituito da una matrice solida, porosa e satura di fluido che scorre attraverso i pori. Il comportamento tempo dipendente è dovuto a due principali meccanismi: i) la viscoelasticità intrinseca delle componenti solide della matrice, ii) la poroelasticità che è un meccanismo estrinseco dovuto al flusso di fluido attraverso i pori. L’obiettivo del presente lavoro è stato lo studio del comportamento meccanico dell’osso corticale con particolare attenzione alla risposta temporale e agli effetti delle lunghezze caratteristiche delle condizioni di prova. La nanoindentazione è una tecnica sperimentale che permette di indagare il comportamento meccanico di un materiale in diverse direzioni. Una prova di nanoindentazione consiste nel sollecitare la superficie del materiale da analizzare con una punta, detta indentatore, costituita da materiale molto duro, ad esempio diamante. I dati di interesse registrati durante una prova sono, in ogni istante, il carico P e l'affondamento h della punta nel materiale, che costruiscono la curva P-h. La prova presenta due fasi distinte di carico e scarico durante le quali viene raggiunto il massimo affondamento o il massimo carico prestabiliti prima di riportare la punta nella posizione originale. A partire da questa curva è possibile estrarre alcuni parametri legati direttamente alle proprietà meccaniche del materiale come il modulo di indentazione e la durezza. Il vantaggio maggiore di questa tecnica risiede nel fatto che, pur operando a scale dell'ordine dei micrometri o dei nanometri e con carichi inferiori ai 300mN, non necessita di requisiti particolarmente stringenti in termini di forma e dimensione dei provini da studiare, risultando quindi uno strumento particolarmente comodo in campo biomeccanico per misurare le proprietà meccaniche dei tessuti biologici a livello microstrutturale e submicrostrutturale. Mentre la viscoelasticità è stata ampiamente studiata e sono disponibili in letteratura diversi modelli analitici per materiali puramente viscoelastici, invece l’interpretazione della risposta meccanica di materiali poroelastici in prove d’indentazione rimane una questione ancora aperta. Va sottolineata la mancanza di prove sperimentali nelle quali vengano eseguite prove su tessuto nelle diverse condizioni sulle quali ci si è concentrati in questo lavoro nonché la non disponibilità di simulazioni ad elementi finiti nelle quali possa essere considerato il comportamento poro-elastico dell’osso corticale e la sua dipendenza dai parametri della simulazione. Inoltre, è necessario evidenziare come, seppur siano state eseguite molte prove in cui il campione veniva mantenuto in condizioni simili a quelle fisiologiche, non siano reperibili esperimenti effettuati con l’ausilio di punte sferiche in luogo della punta Berkovich la quale non è in grado, come qui si vuole dimostrare, di poter identificare il ruolo fondamentale del fenomeno poro-elastico nella determinazione delle proprietà meccaniche del tessuto. Nel presente lavoro viene analizzata la risposta temporale dell’osso corticale in prove di rilassamento e di creep di nanoindentazione eseguite con punte coniche (Berkovich) e sferiche di varie dimensioni. Questo approccio rappresenta una soluzione efficace perché permette, oltre che misurare le proprietà locali e le eterogeneità del materiale, di determinare modalità di deformazione diverse al variare della geometria della punta di indentazione, condizioni di carico e durata della prova. Il presente lavoro è volto a simulare prove di nanoindentazione mediante un modello numerico a elementi finiti sviluppato tramite il codice commerciale ABAQUS che colga gli aspetti di poro-elasticità delle caratteristiche meccaniche dell'osso corticale e che riproduca in maniera adeguata la morfologia della curva P-h con i parametri che la descrivono. A tal fine, viene presentato un confronto tra i risultati delle simulazioni numeriche e le prove sperimentali effettuate da Carnelli et al. (2011), Wang et al.(2006), Bembey et al. (2006), sia per quanto riguarda la corrispondenza tra le curve sia per i valori numerici dei parametri d'interesse. 2. Materiali e metodi La poro-elasticità del tessuto osseo è stata studiata simulando le prove di nanoindentazione in direzioni diverse. In particolare, per il comportamento meccanico del materiale, è stata assunta l'ipotesi di isotropia trasversa secondo la quale, considerando un sistema di riferimento avente asse verticale coincidente con l'asse principale di un osso lungo, il piano radiale-circonferenziale viene inteso come un piano di isotropia. In accordo con questo tipo di simmetria del materiale, le prove di nanoindentazione lungo la direzione assiale sono state simulate mediante un modello agli elementi finiti assial-simmetrico, mentre le prove in direzione trasversale hanno richiesto l'utilizzo di un modello tridimensionale. I dati di letteratura mettono in luce delle marcate differenze tra i valori di sforzo di snervamento a trazione e compressione per il tessuto osseo, rendendo così giustificabile l'utilizzo di un criterio di snervamento pressione-dipendente di tipo Drucker-Prager. Per le simulazioni numeriche sono state sviluppate sia delle prove in controllo di spostamento, per la verifica dei dati in termini di moduli di indentazione e per l’identificazione della geometria della punta ottimale in questo caso, che in controllo di forza, per l’individuazione della dimensione della punta adeguata all’analisi del fenomeno poro-elastico e della permeabilità dell’osso alla scala geometrica considerata. Le prove sono state condotte considerando gli equivalenti computazionali, ideale e arrotondata, della punta Berkovich e punte sferiche di raggio 5,7,10,15,21 e 30μm. Le simulazioni sono state condotte con durata dello step di carico costante pari a 0.1s, durata dello step di holding uguale a 10s e scarico di 1 e 10s. Per le prove di rilassamento i moduli di indentazione drenato e non drenato sono stati ricavati tramite interpolazione dei dati appartenenti al tratto di scarico dal 90% al 40% con la funzione F=a(h-hf)^m proposta da Oliver and Pharr (1992). Il comportamento tempo dipendente della cartilagine viene analizzato tramite interpolazione delle curve di creep con la soluzione analitica per un problema di indentazione di un materiale poroelastico formulata da Agbezuge and Deresiewicz (1974): H^*=A-A/[1+(T^*/T_0 )^P ] Dove spostamento e tempo sono stati normalizzati opportunamente: T^*=√ct/(a(t)) H^*=(h(t)-h(0))/(h(∞)-h(0)) dove c è un parametro costante che dipende dalla permeabilità e dal modulo di taglio. Per quanto riguarda le proprietà elastiche del materiale, le costanti elastiche implementate sono state ricavate dalle prove sperimentali di Carnelli et al.(2011); data l'anisotropia del materiale, i valori analitici di riferimento sono stati poi ricavati introducendo la formulazione esplicita di Delafargue and Ulm (2004). 3. Risultati Dalle prove in controllo di spostamento con l’ausilio della punta ideale e della punta arrotondata è emerso che i moduli di indentazione calcolati nel caso in cui il campione venga messo in condizioni drenate (in cui la risposta è determinata dalla matrice solida) e non drenate (in cui la risposta è dovuta alla componente fluida) non differiscono in maniera significativa. Lo stesso non è stato rilevato, come volevasi dimostrare, per le prove con le punte sferiche in cui i moduli di indentazione riportano una differenza sostanziale tra le due condizioni di prova (Tab 1): Punta M (GPa) Drenato Non drenato Ideale 27.5 29 Arrotondata 25.5 26 R=5μm 23.6 29 R=15μm 23.2 29 R=21μm 23.17 28.7 R=30μm 24 28.5 Tab 1: Moduli di indentazione calcolati a seguito delle simulazioni nelle due condizioni in cui si è posto il campione sottoposto a simulazioni con le versioni computazionali ideale e arrotondata della punta Berkovich e con le punte sferiche qui considerate. Le simulazioni tridimensionali atte al calcolo del rapporto di anisotropia relativamente alle prove eseguite con punte di raggio 15 e 21μm hanno dato luogo a risultati compatibili con i risultati delle prove sperimentali riportati da Carnelli et. al(2011) (Tab.2): Raggio della punta (µm) MT (GPa) MA (GPa) MT/MA HT (GPa) HA (GPa) HT/HA 15 19 23.6 0.8 0.51 0.57 0.89 21 18 23.2 0.78 0.46 0.47 0.97 Tab 2: Moduli di indentazione in direzione trasversale MT e assiale MA, rapporto di anisotropia MT/MA, rigidezza in direzione trasversale HT e assiale HA e rapporto di anisotropia correlato HT/HA, per entrambe le punte considerate in questo caso. Dall’analisi delle curve di creep con punte sferiche è emerso il comportamento poroelastico dell’osso corticale. A differenza di un materiale puramente viscoelastico le curve di creep non si sovrappongono perfettamente per diverse dimensioni della punta di indentazione. La funzione di fitting usata è in accordo con gli andamenti delle simulazioni eseguite per cui è verificata l’applicabilità della soluzione analitica, per un materiale poroelastico sottoposto a indentazione, formulata da Agbezuge and Deresiewicz per i dati relativi all’intera fase di creep come già riscontrato da Oyen et al.(2007) per i dati sperimentali delle prove di nanoindentazione sferica di Bembey et al.(2006). Dal grafico successivo (Fig.2) è possibile verificare l’andamento della costante di tempo T0 associata al fitting rispetto al raggio di contatto associato alle diverse punte. Da questo si può concludere che, per effettuare un’analisi del comportamento poro-elastico del tessuto e una valutazione della permeabilità a questa scala geometrica è opportuno scegliere il raggio della punta nell’intervallo 1-2. Giustappunto, il parametro permeabilità è stato studiato mediante l’uso di una punta sferica di raggio 7μm. Le curve riportate in figura 3 mostrano che gli andamenti relativi alla diminuzione della permeabilità di un decimo e di un ventesimo rispetto al valore iniziale sono piuttosto simili e ciò è indice del fatto che un’ulteriore diminuzione non provocherebbe alcun cambiamento degno di nota nei risultati. Queste considerazioni sono in accordo con i dati sperimentali e teorici indicati in sezione 2.4 dai quali è possibile concludere che un valore di permeabilità inferiore a O(10-25) non è verosimile per l’osso corticale. Osservando infine l’andamento della curva conseguente all’introduzione di un valore di dieci volte superiore rispetto a quello di riferimento in cui il fenomeno poroelastico è eccessivamente veloce, si può affermare che, i valori di permeabilità intrinseca dell’osso corticale ragionevolmente ammissibili oscillano tra O(10-23) e O(10-25) m2. 4. Conclusioni Dai risultati ottenuti si è concluso che se la poroelasticità è un meccanismo fondamentale nella definizione del comportamento tempo dipendente dell’osso corticale, la dipendenza dai parametri costitutivi viene risolta riportando i risultati su una scala temporale normalizzata. Infatti la soluzione analitica per un problema di indentazione di un materiale poroelastico mostra come la risposta temporale dell’affondamento della punta di indentazione sia dipendente dal rapporto T^*=√ct/(a(t)). Ulteriori dipendenze dai fattori dimensionali coinvolgono altri meccanismi dissipativi quali la viscoelasticità intrinseca della matrice solida di cui è composto il tessuto osseo. In questo lavoro di tesi è stata introdotta una nuova metodologia interpretativa del comportamento poroelastico dell’osso corticale, grazie alla versatilità della prova meccanica di nanoindentazione sferica, che consente di effettuare prove caratterizzate da dimensioni geometriche differenti. Tale metodologia di prova, permette quindi, da un lato di rilevare e discernere il comportamento viscoelastico e poroelastico di un tessuto biologico complesso come quello osseo, e dall'altro, di determinare, grazie all'ausilio di una modellazione teorica e numerica i valori numerici di alcuni parametri costitutivi caratteristici di ciascuna scala geometrica. In questo lavoro ci si è concentrati sulla permeabilità, sul modulo d’indentazione e sul rapporto di anisotropia. Quest’ultimo, poiché è stata dimostrata la differenza che intercorre tra le proprietà calcolate tenendo presente la poro-elasticità del tessuto e quelle che invece non comprendevano l’analisi di questo comportamento, deve essere analizzato attraverso l’inserimento di nuove prove tridimensionali che mettano in luce il comportamento del campione anche in condizioni non drenate. L’utilità di questo lavoro si basa quindi, sull’introduzione di un modello ad elementi finiti in cui possa essere enfatizzato il ruolo della componente poro-elastica nella risposta meccanica tempo-dipendente dell’osso corticale, sull’identificazione della geometria e soprattutto della dimensione dell’indentatore che meglio possa interpretare il fenomeno poroelastico di interesse e infine sull’istituzione di una linea guida per l’impostazione di esperimenti atti non solo all’individuazione ma anche al discernimento delle componenti che definiscono il comportamento tempo-dipendente del tessuto.