Home

sviluppo, valutazione e applicazione di metodi

1. 62 92 4 3 Interpolazione locale con funzioni polinomiali L analisi delle 4 13 evidenzia come per x 0 le componenti dello spostamento non siano affatto nulle variando invece in modo parabolica con la coordinata y La definizione della trave di Timoshenko come di una trave incastrata un ap prossimazione lessicale giustificata dall esiguita relativa degli spostamenti so lamente nel contesto della pratica tecnico ingegneristica Ai fini della valutazio ne di precisione degli spostamenti predetti al contrario le condizioni di vincolo devono essere introdotte coerentemente con le ipotesi del problema condizioni sotto le quali le 4 13 valgono esattamente 63 Per riprodurre il problema so no stati dunque imposti i valori di spostamento previsti dalla soluzione esatta su tutti i punti del contorno andando ad indagare la corrispondenza tra soluzione numerica e soluzione teorica nei punti interni I parametri descrittivi del problema trattato nel test sono riassunti in tabella 4 2 il test stato condotto adottando una costellazione costituita da 4 nodi satellite per inodi di interni da tre satelliti per io nodi di bordo e da due per i nodi di ver tice singola cella come metrica di precisione stato considerato il valore qua dratico medio del modulo della differenza vettoriale tra spostamento numerico e spostamento teorico in ciascuno dei punti non vincolati Se Gea si Uist Vi
2. 3 1 1 3 Gradiente Nella generica cella c il vettore gradiente g definito dalle componenti OT 010070 2 0 0 T 0010 O0 0 wW 0 CT 3 24 T 0001077 0 0 2 Le componenti del gradiente sono Gy O T OrE ON db oT Gp ST Oye oyn o l T 3 25 gz 0 T O on o OT dall equazione 3 4 si ha Ie 7 6 a u 0 1007 0 2 0 0 Geno I8 e wy 0 W 20 2 0 0 W OE Ty 8 26 GEM y o ae ODE peo 0 0 X dal punto di vista computazionale vi poi un ulteriore vantaggio nell adozione delle coordi nate locali affini quando i coefficienti della matrice dipendono dalle coordinate essi assumono valori che decrescono al diminuire delle dimensioni della cella e questo peggiora la precisione con cui gli algoritmi disponibili calcolano i coefficienti della matrice inversa 52 3 1 Problemi scalari tridimensionali 3 1 1 4 Equazione costitutiva Il legame tra il vettore flusso di calore q e il gradiente della temperatura g det to equazione costitutiva In un mezzo anisotropo questo legame viene espresso nella forma GE Rae on Kaz Ix Y 2 dyj D 4 he kyy kyz ivo 3 27 qe v y 2 J kzx Key Keel g2 2 y 2 Indicando con K la matrice 373 detta matrice di conduttivit termica della cella c il legame costitutivo pu essere espresso nella forma matriciale de K 8c 3 28 Se il materiale di cui si compone la cella c dal punto di vista della conduzione termica is
3. A 0 0 05 0 1 0 15 0 2 0 2 0 15 0 1 0 05 A mm 0 Figura 6 14 Distribuzione delle differenze nel modulo degli spostamenti a sini stra riportata la distribuzione cumulata dello scarto percentuale riferito al valore UMCM es del FEM uca peau a destra riportata la distribuzione delle differenze UFEM con segno umcm urem 171 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle 172 Conclusioni I modelli numerici tridimensionali subject specific costruiti a partire da dataset di immagini medicali possono costituire uno strumento accurato di predizione del rischio di frattura dei segmenti ossei Ad oggi stata affinata una metodologia di modellazione basata sul Metodo degli Elementi Finiti FEM la procedura ri chiede la precisa definizione topologica della superficie di contorno del segmento osseo e con un processo di modellazione complesso composto da differenti fasi scontornamento descrizione matematica delle superfici di confine generazione della mesh e mappatura delle propriet materiali non tutte completamente au tomatizzate Bench questa procedura abbia dimostrato di poter ottenere un elevata accura tezza nella predizione delle deformazioni di segmenti ossei in vitro la sua ap plicabilit in ambito clinico risulta ancora limitata a causa della sua complessit la generazione di un modello completo a partire dai dati delle
4. 12i pedice c si rende necessario riferendosi all interpolazione locale della costellazione C I3si supposto qui per semplicit che il mezzo abbia comportamento isotropo 97 4 L approccio meshless con il Metodo delle Celle Nell intorno del polo P il calore che attraversa una superficie orientata S di versore normale n evidentemente esprimibile come Qs fancas f ane amy d5 4 19 S S Definita una regione tributaria Qp attorno al nodo polo delimitata dalla super ficie chiusa 00 p il bilancio di nodo impone l equivalenza tra il flusso netto di calore uscente attraverso la superficie Qan e il calore generato Qo all interno del volume da essa racchiuso Qas Qaar 4 20 essendo Qoap r Gx Nz Jy Ny dS Qop 4 21 aQp hy S4 e P x y P R 9 Figura 4 13 Regione tributaria circolare centrata nel nodo polo evi denziato il sistema di riferimento polare locale e la sua relazione con il riferimento cartesiano locale Adottando come regione tributaria una opportuna figura piana elementa re possibile pervenire ad una espressione analitica del flusso di calore netto attraverso la superficie di contorno 98 4 3 Interpolazione locale con funzioni polinomiali Si consideri ad esempio come regione tributaria del nodo polo un cerchio di rag gio R centrato nel polo Il valore di R di principio arbitrario e pu essere conve nientemente assunto pari al val
5. Con il tentativo di superare questi limiti a partire dalla fine degli anni 70 inizi a svilupparsi un approccio differente che costruendo le interpolazioni basandosi unicamente sui nodi non rendesse necessaria la generazione di elementi questa nuova categoria di metodi assunse il nome di meshless talvolta denominati anche meshfree Il primo metodo meshless sviluppato ha nome Smoothed Particle Hydrodynamics SPH sviluppato nel 1977 in ambito astrofisico da Lucy 32 Monaghan e Gin gold 33 fu poi applicato a problemi di fluidodinamica e di meccanica dei solidi problemi di impatto Dall epoca numerosi furono i metodi sviluppati ascri vibili alla categoria dei metodi meshless citando i pi noti e probabilmente tralasciandone in questa arbitraria selezione alcuni e Smoothed Particle Hydrodynamics SPH 32 33 e Non Structured Finite Differences NSFD 34 e Diffuse Element Method DEM 35 e Element Free Galerkin Method EFG o EFGM 36 e Reproducing Kernel Particle method RKPM 37 e partition of Unit PUFEM 38 e Finite Point Method FPM 39 e hp Clouds 40 e Natural Element Method NEM 41 e Meshless Local Petrov Galerkin MLPG 42 e Generalized Finite Element Method GFEM 43 74 4 1 La filosofia dei metodi meshless e Meshless Finite Element Method MFEM 44 e Extended Finite Element Method XFEM 45 Tutti questi metodi condividono il fatto di discretizzare il do
6. di base radiale r ga 1 1 op 1 74 de eon 4 dp V2 r2 4 85 no satelliti 5 lB d dc bp 1 456290 0 913340 0 075073 0 011741 0 329283 0 001377 0 000062 0 000048 0 299443 0 002011 0 000075 0 000025 0 161390 0 001490 0 000068 0 000022 0 079458 0 000898 0 000065 0 000022 0 000139 0 000142 0 000050 0 000031 CON DD O1 amp W Tabella 4 7 Interpolazione con RPF scarto quadratico medio per il pas so griglia pari a 1 4 per le varie funzioni di base radiale in dipendenza del numero dei satelliti della costellazione Un test preliminare ha riguardato la sensibilita della soluzione rispetto al nu mero di nodi utilizzati nella costellazione Assunta una distribuzione di punti regolare con passo di riferimento pari a 1 4 del lato si valutato il parametro di precisione 6 al variare del numero di nodi della costellazione tra 4 e 9 per ciascu na delle funzioni adottate I risultati esposti in tabella 4 7 e rappresentati in figura 4 22 evidenziano come per le funzioni di base radiale B C e D ci sia una sostanziale stabilit rispetto alla numerosit della costellazione adottata per la funzione di base radiale A al contrario per ottenere accuratezze paragonabili a 121 4 L approccio meshless con il Metodo delle Celle 1 600 5 1 400 5 1 200 5 1 000 0 800 5 0 600 5 0 400 5 0 200 5 gt lt Ho e n 0 000 2 wad Sa 5 i 6 num
7. Figura 4 4 Creazione della regione tributaria del nodo polo a sinistra la regione individuata dalla divisione di ogni cella appoggiata al baricentro e ai punti medi dei lati come usualmente effettuato nel CM tradizionale nella figura di destra la regione delimitata in ogni cella da una superficie parallela al lato di cella opposto al nodo polo e passante per il baricentro di cella e cella per cella si calcolano i flussi di calore in transito attraverso la super ficie di confine tra la cella duale del nodo polo e il resto della cella queste grandezze saranno funzione dei valori di temperatura nodale dei vertici della cella esaminata Con riferimento alle notazioni di figura 4 5 il flusso attraverso l elemento di superficie piana orientata A esprimibile come Th 1 Anz Ai A Cr Ax el Ax A K x ix ja T 4 2 dla AT a T 42 J essendo t lo spessore del sistema piano Ac l area della cella c A e A le componenti del vettore area A K la matrice di conducibilit termica della cella c 84 4 2 Meshless a celle locali con il CM Figura 4 5 Flusso attraverso una superficie orientata A nella cella c orien tazione dei vettori area a puro titolo di esempio stata assunta l esatta geo metria della cella 1 di vertici P 5 52 Anz AnyAix i vettori area associati ai lati della cella c di cui alla figura 4 5 e si scrive il bilancio di calore al nodo polo imponendo l equivalenza
8. In letteratura vengono de scritte due strategie Un primo approccio quello di includere nella formulazione degli elementi le propriet meccaniche variabili nello spazio Questo l approccio pi generale che bene si accorda con la sostanziale variabilit delle propriet dei tessuti entro ogni elemento di volume Questo approccio non si pu adottare utilizzando un codice FEM commerciale general purpose richiedendo una formulazione ad hoc e dunque una manipolazione del solutore Inoltre l esigenza di una tale complessit di analisi non stata sinora dimostrata ed ha forse costituito la ragio ne per cui questo approccio risulta essere stato sinora utilizzato in un solo studio 17 Tutti gli altri studi sono basati sull implicita assunzione che le variazioni dei valori CT all interno del volume di ogni elemento siano trascurabili Pertanto le propriet meccaniche del tessuto vengono assegnate agli elementi facendo la media dei valori del campo scalare all interno di ogni elemento Le differenze risiedono nel modo di calcolare questa media all interno degli elementi L approccio pi semplice proposto in letteratura la ricerca per ogni nodo di ciascun elemento del valore associato al pi vicino punto del dataset CT all ele mento viene allora assegnata la media pesata dei valori nodali 18 Una variante proposta in letteratura assegna ad ogni elemento il valore ottenuto mediando le 19 1 La creazione di mode
9. Qoap 27R7tk 0 0 0 1 1 Mg Te 4 35 Il vettore kp 2rR tk 0 0 0 1 1 M Te 4 36 14i caso in oggetto considera per pura compattezza di notazione e senza per questo perdere di generalit una costellazione composta da 4 nodi satellite in questo caso la presenza di un materiale anisotropo porterebbe ad un vettore flusso unitario espresso dalla 4 32 in cui risulte rebbero assenti i termini kzy e kyz se come illustrato nei test di seguito riportati la costellazione a 4 nodi registra buona accuratezza in presenza di un mezzo isotropo per trattare correttamente l anisotropia nella conduzione del calore si rende quindi necessaria una costellazione pi ricca di satelliti 101 4 L approccio meshless con il Metodo delle Celle contiene i valori da inserire nella matrice fondamentale K del sistema risolvente KT 0Q 4 37 alla riga corrispondente all identificativo del nodo polo nella numerazione globa le Le colonne in cui gli elementi di kp vanno collocati corrispondono alle colon ne di indice pari al numero identificativo di ogni singolo nodo della costellazione nella numerazione globale 4 3 1 2 Instabilit Oltre a richiedere una costellazione locale composta da un ben determinato nu mero di nodi satellite l interpolazione polinomiale locale del campo incognito risulta essere poco robusta rispetto alla disposizione geometrica della costellazio ne di punti su cui si appoggia Particolari disposi
10. estesa al problema dello studio delle deformazioni dei corpi elastici in regime di piccoli spostamenti L illustrazione dell implementazione verr completata con un test di accuratez za nel cogliere il campo di spostamenti in un problema di controllo del quale si dispone di una soluzione teorica in forma chiusa i risultati ottenuti vengono comparati con gli analoghi risultati ottenuti nella soluzione del medesimo proble ma con la stessa discertizzazione geometrica e lo stesso ordine di interpolazione da un codice commerciale basato sul Metodo degli Elementi Finiti 58 3 2 Problemi di elasticita lineare tridimensionali 74 8 a T2 94 y 3 9x 47_ aan 10 4 E go 3 ed o 11 4 eae gn O p an SAM eee 12 4 p e TI i I i y 3 8x 5 9 134 pas ip _ N _ W WS e 8 77 ee Le de Le Ne dll 15 ge T T T T T T 2 4 2 2 2 1 8 1 6 1 4 1 2 1 log d Figura 3 6 Analisi dell ordine di convergenza del FEM grafico in scala logarit mica dello scarto quadratico medio e rispetto alla dimensione caratteristica di cella d per le funzioni di test lato cella errore medio FEM T1 T2 0 3333 4 16E 005 1 28E 004 0 2500 1 40E 005 4 25E 005 0 2000 5 89E 006 1 43E 005 0 1667 2 72E 006 7 89E 006 0 1429 1 51E 006 4 34E 006 0 1250 1 07E 006 2 99E 006 0 1111 5 94E 007 1 50E 006 0 1000 4 20E 007 1 16E 006 Tabella 3 2 Valori dello scarto quadratico medio tra soluzione numerica prodo
11. la di appartenenza c tra il vertice n e il nodi di cella adiacente j ciascuno di questi flussi sar espresso in relazione alle temperature nodali della cel la in cui l elemento di superficie in oggetto contenuto Conclusa questa operazione per tutti gli elementi di superficie del poliedro duale possibile scrivere l equazione di bilancio per il nodo in esame yD ee 3 32 ceI n ij essendo il calore generato all interno della cella duale associata ad n analizzando una cella primale alla volta calcolare i flussi di calore attraver so le superfici duali contenute nella cella Componendo le equazioni 3 26 3 27 e 3 31 vengono determinati i termini di flusso delle porzioni di celle tributarie associate ai nodi di cella Questi flussi sono scritti in funzione dei valori delle temperature della cella esaminata Con entrambe le strade si arriva alla scrittura del sistema lineare che defini sce la soluzione del problema le temperature nodali con la prima procedura la matrice fondamentale del sistema risolvente viene scritta riga a riga calcolando gli elementi non nulli corrispondenti ai nodi delle celle che condividono il nodo n in esame i flussi di calore che lasciano il nodo n nella cella c sono espressi in termini dei valori delle temperature associate a tutti i nodi della cella e vengono 55 3 Formulazione quadratica con il Metodo delle Celle considerati i soli flussi che attraversano le superfici di co
12. sin z e sin z 3 34 To z y z e cosx cosy sinx siny 3 35 Come dominio di studio stato considerato un cubo di lato unitario centra to nell origine del sistema di riferimento cartesiano e con le facce ortogonali agli assi del riferimento Il dominio stato suddiviso in tetraedri a 10 nodi creando dapprima una suddivisione in cubi ognuno dei quali stato suddiviso in modo regolare in tetraedri Con questa strategia il lato dei cubi di appoggio ottenuti con la prima suddivisione del dominio fornisce una dimensione lineare caratteri stica associata alla mesh finale di tetraedri Come parametro di errore stato adottato il valore dello scarto quadratico medio tra valori della soluzione analitica e valori della soluzione numerica ai no di della mesh 56 3 1 Problemi scalari tridimensionali Figura 3 4 Esempio di mesh regolare di tetraedri a 10 nodi ottenuta a partire da una suddivisione regolare della geometria La tabella 3 1 riporta i valori dello scarto quadratico medio ottenuti con le due funzioni di test per le varie dimensioni caratteristiche del complesso di celle utilizzato loto di cella Gli stessi valori sono illustrati in figura 3 5 errore medio CM lato cella TI T 0 3333 5 60E 005 1 57E 004 0 2500 1 63E 005 5 15E 005 0 2000 7 99E 006 1 48E 005 0 1667 3 11E 006 7 92E 006 0 1429 1 74E 006 4 36E 006 0 1250 1 40E 006 3 80E 006 0 1111 8 42E 007 1 56E
13. 28 E M Lochmuller O Groll V Kuhn and F Eckstein Mechanical strength of the proximal femur as predicted from geometric and densitometric bone properties at the lower limb versus the distal radius Bone 30 1 207 16 2002 29 M E Taylor K E Tanner M A Freeman and A L Yettram Stress and strain distribution within the intact femur compression or bending Med Eng Phys 18 2 122 31 1996 30 Enzo Tonti A direct discrete formulation of field laws The cell method CMES 2 2 237 258 2001 31 Francesco Cesari Calcolo matriciale delle strutture volume 2 Pitagora Editrice Bologna 1997 32 L B Lucy A numerical approach to the testing of the fission hypothesis Astronomical Journal 82 1013 1024 December 1977 33 J J Monaghan Smoothed particle hydrodynamics Annual Review of Astronomy and Astrophysics 30 1 543 574 1992 34 Q Lin and J Rokne Construction and analysis of meshless finite difference methods Computational Mechanics 37 3 232 248 February 2006 35 Villon P Nayroles B Touzot G Generalizing the finite element method Dif fuse approximation and diffuse elements JournalComputational Mechanics 10 5 307 318 September 1992 36 T Belytschko Y Y Lu and L Gu Element free galerkin methods International Journal for Numerical Methods in Engineering 37 2 229 256 1994 37 W K Liu and Y Chen Reproducing kernel particle methods International Journal for Numerical Methods
14. 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle Figura 5 3 Mediazione delle forze sulla superficie di contorno della ragione tri butaria del nodo P la forza agente sul lato AB assunta pari alla media tra Fi calcolata relativamente all interpolazione degli spostamenti nodali della cel la PQS e F calcolata relativamente all interpolazione degli spostamenti nodali della cella PQR la cella PQR appartiene alla mesh locale del nodo Q Figura 5 4 Regioni tributarie di due nodi vicini distanziate per chiarezza di rappresentazione la superficie di confine tra i due nodi adiacenti stata divisa nelle quattro porzioni di pertinenza di altrettanti tetraedri dei due complessi locali affacciati 136 5 2 Il test esplorativo Nel caso tridimensionale e per disposizioni regolari di punti caso corrispon dente alla situazione del dataset CT una scelta immediata quella di considera re come regione tributaria di ciascun nodo il corrispondente poliedro di Voronoi che in questo caso un parallelepipedo le cui facce sono passanti per i punti medi dei segmenti congiungenti ciascun nodo con i nodi primi vicini secondo le direzioni di allineamento della griglia Costruita la mesh locale di tetraedri la mediazione delle forze avviene sulla superficie di interfaccia tra un nodo e il nodo vicino Questa superficie rientra nella pertinenza di quattro differenti tetr
15. 89 93 pe 12 66 Max Err 272 65 pe 38 37 ue exf 600 400 200 200 400 600 200 400 adduzione 24 y 1 02 x 1 49 R 0 92 RMSE 100 75 pe 11 62 Max Err 263 40 pe 30 38 200 400 600 800 200 400 600 800 flessione 18 y 0 87 x 5 39 R 0 94 RMSE 67 09 pe 11 36 Max Err 159 22 pe 26 96 Figura 6 11 Valori principali di deformazione predetti e misurati nelle sin gole configurazioni di carico 6 4 Comparazione sistematica con modello FEM di riferimento 400 400 300 4 200 4 G 100 4 w w a i T o T T T 1 T T T T T 1 200 100 K 100 200 300 400 300 200 100 100 200 300 400 100 aa PN y 1 045 x 4 639 a y 0 87 x 19 00 R 0 998101 Re 0 990949 200 4 200 4 uz exp uE exp MH 200 y 1 28164x 18 4199 R 0 996406 T T T T T T T 600 500 400 300 200 100 100 200 200 4 T T T T T 1 600 400 200 200 400 600 800 di 400 4 L5 y 0 783 x 3 037 4 Si R 0 99888 pog 600 800 ue exp HE exp Figura 6 12 Valori principali di deformazione predetti e misurati nelle varie configurazioni di carico in alcuni singoli canali di acquisizione quindi quale sia il margine di precisione appurato nella validazione della mo dellazione FEM assunta a riferimento Oggetto del confronto sono i valori principa
16. EnC u cr J3 a ee eae 3 50 che ricordando la relazione 3 41 autorizza la scrittura _ OB End r oB ENC gr OB n6 7 pe die e M u Jt 3 51 Questa espressione consente di calcolare i valori degli elementi del tensore di deformazione e in un punto di coordinate note all interno della cellat c tsi noti come solo termine dipendente dalla cella sia proprio la matrice J legata alla specifica disposizione spaziale dei vertici di cella 63 3 Formulazione quadratica con il Metodo delle Celle 3 2 3 Legame sforzi deformazioni Lo stato di tensione in un punto completamente descritto dalla matrice Orr xy Txz Dia Ciy Fiz 3 52 Ozer Ozy Ozz detta tensore degli sforzi La simmetria di questo tensore consente come per il caso del tensore delle deformazioni di poter adottare la notazione vettoriale che ne riporta le sei sole componenti significative o 3 53 Il legame che correla il tensore degli sforzi o a quello delle deformazioni e descri ve come il materiale si deforma per effetto delle sollecitazioni comportamento macroscopico strettamente legato alla costituzione interna del materiale tipo di legami chimici presenti tipo di struttura molecolare Si parla a tale proposito di legame costitutivo Per un materiale che abbia comportamento elastico lineare ha o De 3 54 Adottando per tensioni o e deformazioni e la notazione vettoriale la 3 54 assu me l
17. Meshfree Methods 52 Irina Ionescu James E Guilkey Martin Berzins Robert M Kirby and Jef frey A Weiss Simulation of soft tissue failure using the material point method Journal of Biomechanical Engineering 128 6 917 924 2006 53 Ken C L Wong Linwei Wang Heye Zhang Huafeng Liu and Peng cheng Shi Meshfree implementation of individualized active cardiac dy namics Computerized Medical Imaging and Graphics 34 1 91 103 2010 182 Image Guided Surgical Planning and Therapy 54 Chun Yang Dalin Tang Chun Yuan William Kerwin Fei Liu Gador Can ton Thomas S Hatsukami and Satya Atluri Meshless generalized finite difference method and human carotid atherosclerotic plaque progression si mulation using multi year mri patient tracking data CMES 28 2 95 108 2 2008 55 Ashley Horton Adam Wittek and Karol Miller Subject specific bio mechanical simulation of brain indentation using a meshless method 2007 56 E Pena M A Martinez B Calvo and M Doblar Application of the na tural element method to finite deformation inelastic problems in isotropic and fiber reinforced biological soft tissues Computer Methods in Applied Mechanics and Engineering 197 21 24 1983 1996 2008 57 H Y Wu K M Liew and T Y Ng Meshless method for modeling of human proximal femur treatment of nonconvex boundaries and stress analysis Computational Mechanics 28 5 390 400 May 2002 58 M Doblar E Cuet
18. RMSE 258 163 RMSE 15 6 9 8 Max err 621 699 Max err 35 5 42 3 Tabella 5 1 Parametri della retta di regressione e valori degli indicatori di accuratezza globale per le predizioni MCM e FEM dati aggregai relativi a tutte le condizioni di carico 147 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle Condizione di Carico LC1 LC2 LC3 LC4 LC5 CM FEM CM FEM CM FEM CM FEM CM FEM Piccolo trocantere U mm 2 3 2 7 1 4 1 7 2 3 2 8 4 2 5 2 1 7 2 2 U mm 0 3 0 4 0 3 0 4 0 6 0 6 0 2 0 2 5 1 6 0 U mm 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 4 0 5 U mm 2 3 2 8 14 1 7 2 4 2 9 4 2 5 2 5 3 6 4 Centro d anca U mm 23 3 9 2 1 2 5 3 3 4 0 5 3 6 6 2 5 3 2 U mm 04 0 5 04 05 08 0 9 04 0 3 6 9 8 2 U mm 0 7 0 9 0 5 0 6 0 7 0 9 0 7 1 0 0 7 0 9 IU mm 3 3 4 0 2 2 2 6 3 5 4 2 5 4 6 7 7 4 8 9 Tabella 5 2 valori dello spostamento al piccolo trocantere e al centro d anca predetti dai modelli MCM e FEM nelle 5 condizioni di carico 148 Capitolo 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle Nel capitolo precedente stata illustrata una prima valutazione esplorativa del l applicazione dell approccio meshless con il CM a celle locali MCM nella pre dizione dello stato di deformazione di segmenti di ossa lunghe La valutazione
19. approssimativo e dimensione eccessiva per riprodurre con accuratezza la reale superficie interna dell osso questo fatto primariamente dovuto al la minore accuratezza del processo di segmentazione del canale Insieme a ci l elevata curvatura della superficie rende maggiormente difficoltosa la creazione di elementi di buona qualit nel descrivere la superficie a meno di non fare ricorso ad elementi di dimensione ridotta questa soluzione non si rende per necessaria atteso il modesto contributo che questi elementi forniscono al comportamento globale del modello essendo posizionati in una zona che per le configurazioni di carico in esame risulta scarsamen te sollecitata L insieme di questi fattori fa s che molti di questi elementi abbiano il centroide collocato all esterno del segmento osseo producendo quindi confronti privi di significato e gli elementi prossimi al punto di applicazione del carico a distanza di meno di 15 mm da esso La valutazione comparativa tra MCM e FEM ha considerato la differenza in valo re assoluto tra i valori principali corrispondenti calcolati dai due modelli riferita al valore predetto dal FEM gl gt pas uca FEMI 6 6 FEM 168 6 4 Comparazione sistematica con modello FEM di riferimento el 10000 8000 J 6000 4000 2000 0 0 0 5 1 1 5 2 400 300 200 100 0 100 200 300 400 Figura 6 13 Confronto tra i valori principali d
20. avendo indicato con E vE E 1 v A gt B 4 101 1 v 1 ni di 2 1 v ELI Calcolando i contributi relativi agli altri tre lati del quadrato sommando op portunamente e considerando il contributo di tutte le funzioni che compongo no l approssimazione locale del campo degli spostamenti si perviene all espres sione X 16d 9y 3x2 4d 3 C 3y2 92 4d 3 A Fan 3 g 1 N 4 yeas azz C B wy 16 ui 3a at i dC i 1 2 2 4 2 2 4 102 componenti della forza netta agente sulla superficie di contorno 02 p in equilibrio con le forze di volume raccolte dalla regione tributaria Qp Essendo P l identificativo del polo nella numerazione globale dei nodi i termini delle sommatorie andranno riportati nelle righe 2P 1 per F Qp e 2P per F Qp alle colonne corrispondenti individuate dagli identificativi dei nodi della costellazione 4 4 2 5 Test numerici Al fine di valutarne la precisione numerica l approccio proposto stato applicato alla soluzione di un problema di elasticit piana noto come trave di Timoshen ko Geometria propriet materiali e carico usati nel test numerico sono riportati in 128 4 4 Interpolazione locale con funzioni di base radiale tabella 4 9 per riprodurre con esattezza le condizioni di validit della soluzione teorica su tutti i punti del contorno sono stati imposti i valori di spostamento previsti dall
21. come la carica dell elettrone o la velocit della luce che sono le stesse per qualsiasi proble ma venga posto in esame quanto le costanti materiali entit costanti nel singolo our usando la formulazione integrale il Metodo dei Volumi Finiti comunque basato sulla formulazione differenziale 26 2 1 Le variabili fisiche problema ma che possono assumere valori di volta in volta differenti in ogni caso specifico La seconda grande categoria di grandezze fisiche costituita dalle variabili fisi che si tratta di grandezze il cui valore definisce una singola particolare condi zione del sistema trattato E una categoria di grandezze assai eterogenea tanto dal punto di vista topologico entita geometriche cui la grandezza viene riferita quanto dal punto di vista funzionale alcune di queste grandezze contribui scono in genere a definire il problema altre sono invece oggetto dell indagine condotta con il modello matematico altre ancora sono grandezze che definisco no lo stato energetico del sistema Per inquadrare il ruolo e l importanza delle variabili si possono adottare due cri teri di classificazione il primo basato sull analisi dell attribuzione alle entit geo metriche del sistema il secondo che indaga invece il ruolo svolto dalla grandezza nel descrivere il sistema studiato 2 1 2 Classificazione in base alle entit geometriche di riferimen to La prima classificazione delle variabili avviene in re
22. e parametri di cui al paragrafo 4 2 2 1 si sono ottenuti i risultati riportati in tabella 5 5 e nel grafico di figura 5 6 che testimoniano la validit in termini di precisione numerica della procedura 5 2 Il test esplorativo L algoritmo MCM stato valutato riproducendo un indagine accoppiata numerico sperimentale condotta con misure di deformazione in campo elastico su un seg mento osseo in vitro 77 un femore umano stato sottoposto a delle prove di carico secondo differenti configurazioni misurando di volta in volta le deforma zioni sulla sua superficie in punti anatomicamente significativi Le prove erano state replicate con un modello FEM costruito con la procedura standardizzata a partire da una acquisizione CT la comparazione tra le misure sperimentali delle 137 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle 2 5 34 lato 3 35 1 9 06E 004 E d 0 5 1 75E 004 I 0 25 4 09E 005 4 51 0 125 1 05E 005 5 a i r i Figura 5 5 Test della trave di Timo q e pa ui shenko per l implementazione me shless a celle locali con simmetrizza Figura 5 6 Test della trave di Ti zione valore dello scarto quadratico moshenko per per l implementa medio 6 per i tra soluzione numerica zione meshless a celle locali con e soluzione teorica nei punti interni simmetrizzazione dati di tabella al dominio al variare del passo della 5 5 in scala lo
23. il calcolo andava a buon fine producendosi stabilmente un risul tato coerente L idea quella di considerare nel bilancio delle forze agenti sul contorno della regione tributaria la media delle forze che si scambiano mutua mente polo e satellite mediando i due valori riferiti alla stessa superficie quando la relazione polo satellite viene rovesciata Facendo riferimento alla figura 5 3 si consideri il nodo P la sua mesh locale di quattro triangoli e la regione tributaria fissando in particolare l attenzione sul lato AB e sia F p la forza agente su AB quando calcolata relativamente alla interpo lazione dei valori di spostamento nodale della cella PQS e sia F la forza agente su AB quando calcolata relativamente alla interpo lazione dei valori di spostamento nodale della cella PQR Nel bilancio delle forze agenti sulla regione tributaria del nodo A relativa mente al lato AB si considera non gi F quanto piuttosto il valore medio 1 2 FE ake ne Si tratta di fatto di considerare lo stesso elemento di super ficie come appartenente ai due campi di interpolazione locale degli spostamenti nodali costituiti dalle celle locali dei due nodi affacciati separati dalla superficie in esame L idea di implementazione immediata nel caso di punti disposti su griglie rego lari generalizzabile con la sola non banale condizione di rispettare il carattere mutuo del rapporto polo satellite tra tutti i nodi vicini 135
24. numerico basato sulla formulazione discreta diretta delle equazioni di campo Il CM scrive le equazioni di campo direttamente nella loro forma algebrica facendo uso di variabili globali riferite agli elementi spaziali e temporali di due comples si di celle tra i quali sussiste una dualit topologica Malgrado sia di ideazione relativamente recente il CM ha conosciuto negli ultimi anni un grande sviluppo venendo applicato in differenti discipline tecnico scientifiche quali l elastostatica l elastodinamica lo studio della trasmissione del calore i problemi di elettroma gnetismo e la dinamica dei fluidi I pregi del Metodo risiedono soprattutto nella semplicit concettuale e nella sua aderenza ella realt fisica e sperimentale la facilit di implementazione e alcuni significativi vantaggi di tipo computazionale rilevati in alcune specifiche applica zioni rendono fondato il considerare il CM una valida alternativa ad altri metodi numerici oggi comunemente utilizzati e nei confronti dei quali ha dimostrato di essere pari e talvolta superiore relativamente a precisione numerica ed ordine di convergenza Il Metodo delle Celle in grado di affrontare sistemi caratterizzati da geome trie complesse e da eterogeneit nelle propriet materiali la formulazione finita che costituisce il cuore del Metodo non trova infatti ostacoli nell affrontare si stemi in cui ogni cella abbia propriet differenti da quelle adiacenti l approccio di
25. te all analisi per canale di acquisizione i valori per alcuni singoli estensimetri sono riportati in figura 6 12 evidenzia da un lato una ricorrente elevata corre lazione tra predizioni e misure sperimentali dall altro una variazione dei coeffi cienti angolari delle rette di regressione Queste osservazioni autorizzano ad affermare che il modello MCM sia in grado di identificare bene la realt sperimentale e che variazioni dei coefficienti angolari delle rette di regressione possano venire attribuite ad un difetto nella definizione delle propriet materiali imputabile a due fattori e gli artefatti dell acquisizione CT la risoluzione finita e non isotropa del dataset produce noti effetti cosiddetti di volume parziale e la non propria assunzione di comportamento isotropo del materiale lega ta all assenza di metodi utilizzabili generali ed accurati per la caratterizza zione dell anisotropia dell osso l accertata anisotropia sia materiale che di struttura viene riprodotta con una disomogeneit spaziale di modulo di elasticit sia nel modello MCM che nel modello FEM 164 6 4 Comparazione sistematica con modello FEM di riferimento 800 gt ue MCM 1000 800 600 400 800 valori aggregati 209 y 0 99 x 4 21 R 0 93 RMSE 78 29 ue 9 03 400 Max Err 272 65 ue 31 45 600 800 Figura 6 10 Femore 2506 rappresentazione aggregata delle deformazioni prin cipali pred
26. termica k del mezzo dc kgc 4 72 Nell intorno del polo P il calore che attraversa una superficie orientata S di versore normale n dunque esprimibile come Qs fan dS iw Ne qy Ny dS 4 73 S S Definita una regione tributaria Qp attorno al nodo polo delimitata dalla super ficie chiusa O00 p il bilancio di nodo impone l equivalenza tra il flusso netto di calore uscente attraverso la superficie Qaq e il calore generato Qo all interno del volume da essa racchiuso Qop Qanp 4 74 essendo Qaop r Gx Mz Jy Ny dS Qop 4 75 aQp 4 4 1 2 Bilancio al nodo polo Attorno al nodo polo consideriamo come regione tributaria un quadrato di lato 2d con baricentro corrispondente al polo stesso il valore di 2d arbitrario e viene qui assunto pari alla met della distanza media dei nodi satellite dal nodo polo Non appare superfluo ricordare che in presenza di sorgenti distribuite occorre che l insieme delle regioni tributarie costruite attorno ai vari punti di volta in volta considerati come poli di una costellazione locale copra l intero dominio Il calcolo dei flussi di calore attraverso ogni lato della regione tributaria risulta cos semplificato Con riferimento alla figura 4 21 si ha infatti essendo t lo spes sore del sistema piano che per il lato AB il vettore area corrispondente risulta si suppone qui per semplicit che il mezzo abbia comportamento isotropo 118 4 4 Interp
27. 0 0 02 eee eee eee 6 2 1 6 2 2 Attribuzione delle propriet meccaniche nell MCM L individuazione dei valori puntuali di deformazione 6 3 Validazione con le misure sperimentali 6 4 Comparazione sistematica con modello FEM di riferimento Conclusioni Bibliografia II IV Introduzione La predizione del rischio di frattura nei tessuti ossei un tema dalle molteplici ricadute nella pratica clinica sia nell ambito dei processi di diagnosi che in quelli di prognosi e di definizione di terapie L accurata predizione del rischio di frat tura pu infatti risultare importante nel trattamento dei pazienti osteoporotici nella pianificazione e nella gestione delle ricostruzioni scheletriche in oncologia pediatrica cos come nella pianificazione preoperatoria degli interventi di chirur gia protesica articolare La stima del rischio di frattura avviene oggi valutando con metodologie radio grafiche bidimensionali il contenuto minerale dell osso in specifici siti anatomici ritenuti significativi il rischio di frattura ottenuto tramite inferenze statistiche rispetto al contenuto minerale ottenendo predizioni il cui grado di attendibilit dell ordine del 70 La frattura per un evento che si sviluppa nell osso a livello di organo e in dipendenza di fattori sia interni quantit di osso sua distri buzione e struttura spaziale che esterni i carichi In considerazione di ci risul
28. 006 0 1000 4 68E 007 1 14E 006 Tabella 3 1 per le due funzioni di test 57 Valori dello scarto quadratico medio tra soluzione numerica prodotta dal CM e soluzione teorica al variare del lato di cella della mesh 3 Formulazione quadratica con il Metodo delle Celle log 7 7 A T2 8 4 y 4 1x 4 4 ur hi 9 4 Pr ait a Pi wae e 77 10 a ge Ka ger ee ger pe 7 11 ae oe we we 4 0 77 12 pe sg a A e T1 13 4 e ral y 3 9x 5 6 a a ga Ss 14 e 07 15 F T T T T T T l 2 4 2 2 2 1 8 1 6 1 4 1 2 1 log d Figura 3 5 Analisi di convergenza rappresentazione in scala logaritmica del lo scarto quadratico medio per le funzioni di test 3 34 in blu e 3 35 in rosso rispetto alla dimensione rappresentativa del lato di cella Gli stessi test utilizzando lo stesso insieme di tetraedri e con le stesse condi zioni al contorno sono stati replicati con il software commerciale Ansys utiliz zante il Metodo degli Elementi Finiti con tetraedri ad interpolazione quadrati ca La tabella 3 2 riporta i valori dello scarto quadratico medio ottenuti con le due funzioni di test per le varie dimensioni caratteristiche del complesso di celle utilizzato loto di cella Gli stessi valori sono illustrati in figura 3 6 3 2 Problemi di elasticit lineare tridimensionali In questo paragrafo l implementazione del Metodo delle Celle con interpolazione quadratica in domini tridimensionali verr
29. 1 4 PI c PJ invertendo la quale si ottiene il vettore dei coefficienti del polinomio interpolante espresso in funzione delle temperature nodali Ah Qi dj ak Al 3 3 3 1 NO NI DI dio 0 si o NODO 0 0 0 1 Dir i DI 0 0 oOooeooo co 0 4 0 O o o o o blIOOIIO lla linearit del legame tra coordinate affini locali e indifferente la scelta del sistema di riferimento 51 0 00 0 00 i 4 00 o 0 o 7 eG ala 3 21 4 4 0 0 0 0 de Ne 4 0 0 uae coordinate cartesiane locali rende 3 Formulazione quadratica con il Metodo delle Celle ovvero acts 3 22 L adozione delle coordinate affini locali offre dunque un significativo vantaggio i coefficienti della matrice dell equazione 3 20 sono costanti indipendenti dalla cella dalla sua geometria e dalla sua disposizione spaziale Adottando il riferi mento cartesiano globale i coefficienti della matrice sarebbero stati invece funzio ne delle coordinate dei nodi della cella e questo avrebbe obbligato al calcolo della matrice inversa per l equazione 3 21 per ogni cella con conseguente maggiore costo computazionale della procedura L espressione 3 19 pu dunque essere scritta nella sintetica forma FEMALE Cp GCE E COL 3 23 che permette di calcolare il valore della temperatura in un punto P interno alla cella c noti i valori delle sue coordinate n nel sistema di riferimento locale
30. 41500 pe num E 1000 y 0 97x 1 50 R 0 95 il Fem y 0 93x 1 88 R 0 94 e cm 1500 500 ME GXp 1000 Figura 6 9 Valori principali di deformazione sperimentali asse x e numerici asse y ottenuti con i modelli MCM e FEM valori aggregati per tutti gli 8 femori e tutte le 6 condizioni di carico ratezza delle predizioni del modello MCM sono risultati paragonabili ai valori caratterizzanti le predizioni dei modelli FEM di riferimento 163 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle strain correlation MCM FEM R 0 94 0 95 pendenza 0 93 0 97 intercetta 1 88 1 50 RMSE 91 16 83 11 RMSE 7 92 7 22 Max err 552 31 478 75 Max err 47 99 41 59 Tabella 6 1 Parametri della retta di regressione e valori degli indicatori di ac curatezza per le predizioni MCM e FEM dati aggreganti i valori relativi a tutti i femori in tutte le condizioni di carico Nei grafici delle figure 6 11 e 6 12 sono riportati a titolo di esempio i con fronti predizioni misure nelle varie condizioni di carico relativamente ad uno dei femori in esame che per il suo carattere di generalit nel risultato dell analisi assume in questo ambito valore paradigmatico come caso di riferimento per le considerazioni sulle prestazioni del modello MCM L analisi delle singole condizioni di carico riportate in figura 6 11 unitamen
31. Di fatto le equazioni costitutive costituiscono un ponte tra le variabili di confi gurazione e le variabili di sorgente di uno stesso campo Il requisito di uniformit del campo pone dei limiti precisi alla validit delle equa zioni costitutive il problema viene affrontato in modo diverso dalla trattazione differenziale e quella discreta L approccio differenziale scrive le equazioni rela tivamente ad intorni infinitesimi del punto situazione in cui qualsiasi funzione 41 2 Il Metodo delle Celle approssimabile al primo ordine e di conseguenza il suo gradiente risulta uni forme La formulazione discreta considera invece porzioni di sistema piccole ma finite all interno delle quali si considera il campo uniforme Questa strada segui ta dal Metodo delle Celle accetta a priori un approssimazione cui si comunque costretti all atto di discretizzare le equazioni ottenute con l approccio differenzia le che vengono forzate non senza accettare approssimazioni ad avere validit su porzioni di sistema di estensione finita 2 6 3 I diagrammi di Tonti Per la formulazione di un modello di analisi di un sistema fisico viepi quando si intenda fornirne una formulazione discreta essenziale individuare le gran dezze coinvolte e gli enti geometrici cui esse sono riferite Secondariamente in termini logici necessario rilevare le relazioni che sussistono tra le grandezze tra di loro e nella loro eventuale evoluz
32. Innanzitutto l interpolazione locale pu essere espressa nella forma T P T zy a br cy dr ey l xy y 4 27 aao 2 Imponendo che la funzione assuma nei nodi della costellazione i valori del campo ad essi associati si ottiene il sistema di equazioni lineari lay zi yi a T 1 2 y2 x y b 13 1 z3 ig at y lt c T 4 28 1 z4 ys ci vi d T 1 5 ys 25 yu le Ts sinteticamente esprimibile in forma matriciale come Me ac Te 4 29 ove con il pedice si inteso indicare il riferimento alla costellazione di punti esaminata 100 4 3 Interpolazione locale con funzioni polinomiali Il vettore dei coefficienti del polinomio interpolante viene quindi ad essere espres so dalla ac Mz Te 4 30 Il vettore gradiente pu venire scritto come _ Gel Gad 0 Ze 0 e n 010 alt wy e il flusso unitario assume quindi la forma Is 26 gut 10 2z 0 qe oe k 01 0 2y ac 4 32 Se il materiale fosse anisotropo alla costante k andrebbe sostituita la matrice krz 4 33 j Kyy Il flusso netto attraverso la superficie di contorno della regione tributaria del polo essendo n n n il versore della normale alla superficie assume la forma Qanp i n qe dS aNp 0 10 2x 0 f al 01 0 A ac dS 25 ee Rcos 6 Rsin 6 i RIO Ses i ac d0 0 4 34 0 0 1 0 2Rsin 0 2rR tk 0 0 0 1 1 a essendo t lo spessore del sistema piano Utilizzando la 4 30 si ottiene immedia tamente
33. Instron Canton MA USA con una piattaforma che consente l inclinazio ne del campione secondo i diversi angoli di prova corrispondenti alle condizioni di carico per evitare la trasmissione di indesiderate componenti orizzontali di carico il collegamento alla cella di carico stato realizzato con slitte a basso at trito Il protocollo di misura ha permesso di quantificare la linearit del compor tamento del materiale i fenomeni viscoleastici e il manifestarsi di deformazioni permanenti secondo i dettagli gi esposti al paragrafo 5 2 1 152 6 1 L esperimento di riferimento 6 1 2 I modelli FEM La realizzazione dei modelli FEM basata sui dataset CT di ognuno dei femori ha seguito una procedura consolidata 77 con la segmentazione di ogni data set software AMIRA 3 1 1 Mercury Computer Systems Inc USA si definito un modello geometrico NURBS di ogni femore Geomagic Studio v 7 Rain drop Geomagic Inc USA con sui stata costruita una mesh non strutturata di tetraedri a 10 nodi generata automaticamente da un algoritmo advancing front HyperMesh v 7 Altair Engineering Inc USA In corrispondenza del centro di ciascuno estensimetri stato posizionato un nodo con una procedura di registra zione con il riferimento sperimentale cfr paragrafo 6 1 3 Ogni modello FEM stato tagliato e vincolato a incastro a livello della base di cemento della prova sperimentale La dimensione media degli elementi superficiali di
34. ambito della ricerca La progettazione di una protesi ad esempio non pu certo prescindere dalla conoscenza delle sollecitazioni a cui sar sottoposta cos come pure la previsione gli stress sopportabili da un arto ricostruito importante ai fini di un appropriata programmazione dell attivit di riabilitazione potendo controllare i rischi di frattura rispetto alle sollecitazioni attese La predizione del rischio di frattura assume poi rilievo in ambito clinico quale strumento diagnostico di prognosi e di cura 1 1 Scopi e ragioni dell analisi La misura delle sollecitazioni meccaniche in vivo ovvero su pazienti in vita non fattibile senza evitare invasivi interventi chirurgici operazione che peraltro comunemente ritenuta eticamente inaccettabile Il solo modo per poter fare una previsione non invasiva degli stress meccanici in vivo la costruzione di un mo dello numerico L estrema variabilit sia in termini di morfologia che di pecu liare distribuzione spaziale delle prorpiet fisiche e meccaniche dei segmenti os sei impongono nella prospettiva delle applicazioni cliniche la manipolazione di modelli specificamente costruiti sul singolo caso da studiare tali modelli sono solitamente detti modelli subject specific L individuazione di un metodo per poter simulare con attendibilit le sollecita zioni prodotte comporta poi ovviamente tutti i vantaggi gi noti per ogni settore ingegneristico ridurre notevolmente le oner
35. arbitraria delle celle primali e costituiscono le regioni tributarie associate ad ogni vertice o in termini generali 0 cella del complesso primale Sotto questo profilo il complesso duale di celle pu essere visto come una struttura di suppor to del sistema primale volta a collocare le grandezze coinvolte nell equazione di bilancio La stessa equazione di bilancio associata ad una cella duale la solu zione di tutte le equazioni di bilancio associate a tutte le celle del sistema duale consente di determinare in ultima analisi le variabili di configurazione associate al sistema primale Merita a questo punto sottolineare alcune propriet immediate delle equazio ne di bilancio e indipendenza rispetto alla forma della regione tributaria di riferimento e indipendenza rispetto all estensione della regione tributaria di riferimento 79 4 L approccio meshless con il Metodo delle Celle e relativamente ad un sistema in equilibrio risulter in equilibrio anche qua lunque sua porzione arbitrariamente scelta questo significa che relativa mente ad un sistema di equilibrio possibile scrivere n equazioni di bi lancio riferite ad altrettante regioni estratte dal sistema originale queste regioni in particolare possono sovrapporsi in tutto o in parte senza che la validit dei bilanci ad esse riferiti ne traggano danno alcuno Queste considerazioni riassuntive sul Metodo delle Celle e sui suoi strumen ti sugge
36. associate ai poli di volta in volta considerati ottenuti dalla suddivisione delle celle della mesh locale in a b e c sono illustrate mesh locali di elementi triango lari ottenute considerando rispettivamente 4 8 e 6 nodi satellite in d illustrata la situazione di un nodo di bordo Esclusa la natura globale della costruzione della mesh e del bilancio al nodo l approccio non presenta significativi elementi di novit le celle del comples so locale costituiscono il supporto per l attribuzione delle grandezze coinvolte nel problema fisico trattato e identiche restano le relazioni che portano alla co struzione del bilancio locale ai nodi di cella Nei due paragrafi seguenti questo approccio meshless a celle locali verr presentato nelle sue applicazioni al pro blema della condizione termica in regime stazionario e al problema dello studio della deformazione di sistemi piani in regime di piccoli spostamenti due ambiti assunti come riferimento in tutte le indagini esplorative di questo elaborato 82 4 2 Meshless a celle locali con il CM 4 2 1 Conduzione termica in regime stazionario Si assuma che in un sistema di geometria propriet fisiche conducibilit termica k note siano stati posizionati un conveniente numero di punti o nodi sia al l interno che sul contorno avendosi in particolare un nodo in corrispondenza di ogni eventuale vertice La posizione dei nodi sia nota rispetto ad una riferimento cartesiano global
37. baricentro G x6 yc zc della superficie A per il vettore area a della superficie stessa Usando la notazione vettoriale per il tensore degli sforzi si pu scrivere dunque T Gy 00 A I 10 a 0 ap a 0 T 0 0 a 0 ay ae rA 3 61 4 G ove si inteso adottare il pedice G per indicare che le componenti del vettore o sono calcolate nel baricentro G della superficie A 65 3 Formulazione quadratica con il Metodo delle Celle 3 2 5 L equazione fondamentale il bilancio delle forze L equazione fondamentale dell elastostatica lega direttamente le forze esterne con gli spostamenti ai nodi delle celle Si ottiene componendo opportunamente le relazioni sopra descritte Pape te Stop aL 3 62 Per ottenere la relazione fra forze e spostamenti ai nodi di una generica cel la c si scrive un semplice bilancio delle forze pertinenti a ciascun nodo Questo bilancio viene scritto relativamente alle celle del complesso duale ovvero a quel le regioni del dominio che la divisione in celle primali associa naturalmente ad ogni nodo queste regioni possono essere viste come le porzioni di dominio di competenza di ogni nodo Le celle duali costituiscono delle regioni tributarie dei nodi del sistema primale poich relativamente ad esse che il bilancio delle forze viene scritto contemplando tutte le forze da esse raccolte Di seguito si descrive il bilancio dalla prospettiva della generica cella primale Nella creazi
38. centro griglia digitalizzata e riportata sul modello di ciascun estensimetro A causa della presenza di elevati gradienti sia di densit che di stress nella maggior parte delle posizioni in cui sono stati posizionati gli estensimetri le tensioni principali usate per il confronto con le misure sperimen tali sono state mediate sui nodi di superficie corrispondenti all area sensibile di ciascun estensimetro 5 2 6 2 Modello MCM Nel modello MCM viene determinato il vettore spostamento in ogni punto del dataset CT All interno di ogni voxel quindi individuabile un campo di sposta menti continuo interpolando con una funzione trilineare i valori di spostamento ai vertici I valori di deformazione nei punti corrispondenti agli estensimetri so no stati quindi calcolati derivando il campo di spostamenti del voxel contenen te il punto corrispondente all estensimetro i valori principali di deformazione adottati sono quelli corrispondenti alle direzioni di tensione minima e massima conformemente all procedura sopra descritta a proposito del modello FEM 5 2 6 3 Metrica di accuratezza L accuratezza delle predizioni numeriche rispetto alle misure di deformazione ricavate sperimentalmente stata verificata nel modo seguente e stata realizzata una regressione lineare tra parametri omologhi numerici vs sperimentali la bont della predizione in questo caso individuata dal valore di correlazione R e dai parametri coefficiente angolar
39. dall esigenza di garantire delle propriet formali prodotte dall orientazio ne interna Gli istanti sono orientati come pozzi ovvero sono considerati positivi gli intervalli che arrivano all istante J e negativi quelli che escono dall istante J 2 3 2 Orientazione esterna L orientazione esterna anche per gli elementi temporali pi complessa Consideriamo due intervalli di tempo consecutivi T e Tt come indicato in fi gura 2 4 e consideriamo un istante interno a ciascun intervallo ad esempio gli 33 2 Il Metodo delle Celle istanti medi indicati rispettivamente con le notazioni T ed I Questi istanti ven gono detti istanti duali poich si pongono come corrispondenti agli intervalli di tempo definiti da coppie di istanti primali L intervallo di tempo T che decorre tra due istanti duali prende il nome di intervallo duale Vi un evidente corri spondenza dualit fra istanti primali e intervalli duali e fra intervalli primali e istanti duali L orientazione interna di un istante primale induce naturalmente un orientazio ne del corrispondente intervallo duale Questa appenda descritta l orientazione esterna di un intervallo che corrisponde per analogia ad una sorta di trazione o compressione dell intervallo Analogamente l orientazione interna di un intervallo primale trasferisce sull i stante duale corrispondente una direzione di attraversamento che ne costituisce l orientazione esterna T int
40. diagramma perde evidentemen te la sua tridimensionalit e si riduce ad una rappresentazione piana dei legami tra le grandezze coinvolte Le figure 2 7 e 2 8 riportano il diagramma di Ton ti per il problema elastico riassumendo con evidente sinteticit le relazioni fon damentali della meccanica dei continui nella formulazione differenziale e finita rispettivamente 43 2 Il Metodo delle Celle configuration variables space inner orientation 1P initial k displacement displacement gradient V h Vit incompatibility tensor D ne ViuF Elastostatics differential formulation Lam Navier GAU A G V V u f cubic dilatation pressure K p Hooke CD source variables space outer orientation I L_a Tij 2G ei strain deviator k hij ajkh aj metric tensor 5 hij hij 3 AV o 1 ij ij Oi Yeg V u y k stress deviator stress tensor iv force for unit volume 28 stress tensor K bulk modulus G shear modulus Ref Fung Y C Foundation of Solid Mechanics Prentice Hall 1965 p 129 Xki stress functions Figura 2 7 Diagramma di Tonti per la formulazione differenziale dell elestostatica 2 6 Impostazione del problema con il Metodo delle Celle Elastostatics algebraic formulation configuration variables small displacement th
41. il Metodo delle Celle 4 4 2 3 Sforzi Per un materiale isotropo a comportamento elastico lineare il legame tra sforzi e deformazioni come noto descritto dalla relazione lineare lov 0 A ee de ae gh Ey 4 95 1 rv 1 v Try 0 O 5 Vay da cui tramite le 4 94 4v y y 1 y yi z 1 E COE 1 n2 a e Cio 1 y y e 2i E 1 7 Pe ot Baw 1 y i 2 2 E a H 1 p a 4v xu 2 1 yu yi xu xi E 4 96 5 1 v o 40 0 l 0 e ma E AER 2 1 v l gt 4 y yi 1 y y 2 2 1 v E i 2 1 v 4 4 2 4 Bilancio delle forze al nodo Non appare superfluo ricordare che in presenza di sorgenti distribuite occorre che l insieme delle regioni tributarie costruite attorno ai vari punti di volta in volta considerati come poli di una costellazione locale copra l intero dominio Come regione tributaria del nodo polo si assuma un quadrato di lato 2d con baricentro corrispondente al polo stesso la dimensione del lato come detto in precedenza assolutamente arbitraria e verr qui assunto convenzionalmente pa ri al valore della distanza media dei nodi satellite del nodo polo 66 va ricordato per che in presenza di sorgenti distribuite forze di volume affinch nessuna loro porzione venga tralasciata nella scrittura dei bilanci l insieme delle regioni tributarie dei vari nodi deve coprire tutto il dominio in esa
42. immagini dia gnostiche richiede infatti un notevole impegno sia in termini di tempo richiesto tipicamente alcuni giorni che di competenze specialistiche continuativamente dedicate ad ogni singolo studio fatto che rende praticamente impossibile una sua applicazione su vasta scala con tempi contenuti Scopo della tesi era l esplorazione di un approccio di modellazione alternati vo in grado di conservare i requisiti di robustezza accuratezza e generalit ma caratterizzato da un maggiore grado di automazione in grado cio di approntare e completare un analisi con il minimo grado possibile di intervento specializza to L indagine si concentrata su implementazioni innovative del Metodo delle Celle CM un metodo numerico basato sulla formulazione discreta diretta delle equazioni dei campi fisici l uso di variabili globali riferite agli elementi spaziali e temporali di due complessi di celle primale e duale porta alla scrittura delle equazioni di campo direttamente nella loro forma algebrica La scelta di con centrare lo studio sul Metodo delle Celle legata a pi fattori da un lato alcuni elementi caratterizzanti il CM semplict concettuale e implementativa attitudi ne ad affrontare sistemi geometricamente complessi e disomogenei assenza di singolarit prodotte da sorgenti o da vincoli concentrati lo rendono adatto alle 173 applicazioni nell ambito della biomeccanica ortopedica dall altro vi una con cordanza con
43. in cui esse ricadono Le forze di volume sono note quando sia noto il volume della regione tributaria e le propriet la massa volumica nel caso della forza peso del materiale di cella 66 3 2 Problemi di elasticita lineare tridimensionali Le forze di superficie sono anch esse note quando sia nota la geometria e l orien tazione della superficie o x y z ovvero o 7 nella cella c in esame da stimare nel baricentro della sueprficie Dalla 3 60 si ha infatti T 0 G Au Tr o G Ar 3 63 Ti 0 G As La notazione c G suggerisce che a sia calcolato nel baricentro della superficie Ap Il bilancio delle forze sulla porzione di cella duale del nodo h contenuta nella cella c viene allora scritto nella forma y 5 E E Ta 3 64 cES h essendo S h l insieme delle celle che insistono sul nodo h ossia quelle che hanno h tra i vertici di cella Le forze T possono essere espresse in funzione del vettore degli spostamenti nodali u componendo opportunamente le espressioni 3 41 3 45 3 61 e 3 55 Scrivendo il bilancio su tutte le porzioni di cella duale della cella primale c si ottiene un espressione del tipo F ku 0 3 65 ossia in modo pi usuale F ku 3 66 dove F riassume i termini delle forze di volume e delle forze applicate diretta mente ai nodi Il sistema 3 66 prende il nome di sistema fondamentale e la matrice k detta matrice di rigidezza locale e per la cell
44. l origine del sistema globale nel nodo polo Ne Ung gt Wine gt ri 4 86 vey gt _ Wine gt ri i 1 essendo r il modulo del vettore posizione del punto P x y dal nodo P della costellazione C n P PI ye z y 4 87 Come funzione di base radiale si adotter a scopo di illustrazione la funzione FESSER 4 88 I pesi delle due interpolazioni RBF vengono espressi in termini componen ti nodali del vettore spostamento che costituiscono l incognita del problema applicando la condizione di interpolazione e ottenendo i sistemi di equazioni lineari Bewe uc 4 89 da cui We D_ Uc 4 90 dove si sono raccolti le componenti di spostamento nodale e i pesi nei vettori Ul Wiu Uy Wly wae wee ae 4 91 UN WNu UN UNv avendosi quindi r11 0 d 712 0 ms Oc 0 0 ria 0 Velie sce 0 ri e 2 i 4 92 prno O rewva O gt Olry w 0 0 b Tna 0 rN 0 Tw 124 4 4 Interpolazione locale con funzioni di base radiale Come per il caso dell interpolazione polinomiale le componenti dello spo stamento verranno opportunamente derivate per esprimere le componenti del tensore delle deformazioni con cui a mezzo della legge costitutiva del materia le verranno espresse le componenti del tensore degli sforzi definita una regione tributaria attorno al nodo polo le componenti della forza risultante agente sulla sua superficie di contorno sar calcolata integrando opportuname
45. misura ottenuta dal rapporto tra dimensione dello spigolo radiale e numero di divisioni una misura lineare significativa della mesh 69 3 Formulazione quadratica con il Metodo delle Celle Hi Figura 3 8 Sfera cava pressurizza Figura 3 9 Sfera cava pressurizza ta geometria in esame ta esempio di mesh ui previsto dalla formula 3 68 Come parametro di precisione della soluzio ne numerica si considerato lo scarto quadratico medio rispetto alla soluzione teorica Ni num __ qt 2 essendo N il numero dei nodi non vincolati La tabella di figura 3 10 riporta i valori di ottenuti con il codice CM e con l analisi FEM accoppiata per vari valori della dimensione lineare caratteristica della mesh l I medesimi valori sono riportati in scala logaritmica nel grafico di figura 3 11 dove si evince oltre alla piena concordanza dei risultati ottenuti dai due metodi e alla loro precisione nel predire i valori teorici di spostamento anche la stretta corrispondenza in termini di ordine di convergenza sussistente tra CM e FEM 70 3 2 Problemi di elasticita lineare tridimensionali lato errore medio 6 CM FEM 1 0 0182709 0 0182706 0 5 0 0113420 0 0115367 0 25 0 0031108 0 0031523 0 125 0 0008425 0 0008515 Figura 3 10 Sfera cava pressurizzata rappresentazione del campo di spostamenti e valori dell errore quadratico me
46. segmenti ossei in vitro La metodologia affinata si basa sulla precisa rappresentazione topologica della superficie dell osso fatto che implica un processo complesso composto da pi fasi l individuazione della morfologia del segmento osseo mediante la defi nizione dei suoi contorni in ciascuna delle immagini del dataset CT operazione denominata segmentazione ad oggi ancora al pi semi automatica la gene razione di una descrizione matematica della superficie di contorno dell osso la generazione della mesh attualmente condotta con un software automatico e la mappatura delle propriet meccaniche all interno della mesh Bench notevol mente progredito negli ultimi anni questo processo accusa ancora un importante limite rappresentato dalla complessit procedurale la generazione di un model lo completo a a partire dai dati delle immagini diagnostiche richiede un notevole impegno sia in termini di tempo richiesto tipicamente alcuni giorni che di com petenze specialistiche continuativamente dedicate ad ogni singolo studio Questo rende praticamente impossibile un applicazione su vasta scala con tempi compa tibili con quelli delle esigenze cliniche Nella letteratura di settore sono stati negli ultimi anni esplorati seppur a livello coacervale approcci alternativi nell intento di superare i limiti costituiti da que sta onerosit di modellizzazione La presente tesi formula e analizza appunto alcune nuove ipotesi metodologich
47. spaziali 2 2 1 Orientazione interna degli elementi spaziali 2 2 2 Orientazione esterna degli elementi spaziali L orientazione degli elementi temporali 2 3 1 Orientazione interna LL 2 3 2 Orientazione esterna LL Attribuzione delle variabili agli elementi geometrici licomiplessindi celle xegsitogu ea e pa ani Impostazione del problema con il Metodo delle Celle 2 6 1 Equazioni di struttura 2 6 2 Equazioni costitutive 02 4 sa secon Ple dadi di i eS 2 6 3 Idiagrammi di LOM Ss ser elio kpin alsin amp paci kon i 27 3 Formulazione quadratica con il Metodo delle Celle 3 1 Problemi scalari tridimensionali 00 3 1 1 Interpolazione quadratica 00048 Dili Testdi conyerge a so Id yak scar il Gos iii 3 2 Problemi di elasticit lineare tridimensionali 3 2 1 Interpolazione quadratica del campo degli spostamenti 3 2 2 Spostamenti deformazioni 1003 54 x Seed ria 3 2 3 Legame sforzi deformazioni onda sl anna 3 2 4 Legame forze sforzi a fon oono a 3 2 5 L equazione fondamentale il bilancio delle forze 3 2 6 Test HUMETICO 0 sa ae a e A e he ke 4 L approccio meshless con il Metodo delle Celle 4 1 La filosofia dei metodi meshless 004 4 1 1 Utilizzo di approcci meshless in ambito biomeccanico 4 2 Meshless a celle locali con il CM 4 2 1 Conduzione termica in regime sta
48. studio della deformazione di un segmento osseo in vitro pro blema di riferimento per la valutazione di accuratezza del modello in un conte sto affine a quello reale sia per complessit geometrica che per comportamento materiale 76 L obiettivo quello di accertare la validit dell implementazione nello speci fico contesto delle applicazioni di meccanica dell osso verificando l accuratezza dei risultati rispetto sia alle misure sperimentali che ai valori predetti da un mo dello FEM accoppiato di riferimento La prospettiva dello studio come detto l applicazione clinica ambito che impone metodologie di indagine e robuste rispetto alle singolarit e agli artefatti della sorgente di informazio ne circa geometria e prorpiet materiali e automatiche nella loro realizzazione riducendo idealmente eliminando le fasi di manipolazione manuale altre dalla definizione di carichi e vincoli e generali nella loro applicazione non contenendo elementi restrittivi l appli cabilit a specifici siti anatomici Questi requisiti indirizzano lo studio applicativo verso una costruzione del mo dello numerico basata direttamente sui dati delle immagini medicali i dataset CT specificamente che sono lo strumento di indagine in grado di fornire infor mazioni circa la morfologia e la distribuzione di densit dei tessuti ossei tratta ti Peraltro uno dei maggiori limiti della procedura oggi di riferimento basata sul FEM rappresentat
49. sviluppo teorico e come altre altre elaborate tecniche siano attingibili dalla letteratura 46 e il calcolo delle derivate risulta per alcuni algoritmi significativamente labo rioso e peri metodi che usano i punti di Gauss per integrare la formulazione for te della forma differenziale il numero di punti richiesto in alcuni casi elevato e le funzioni di forma hanno spesso continuit maggiore di C con riduzione Figura 4 1 Discretizzazione con approccio meshless nodi e domini di influenza Q 7 76 4 1 La filosofia dei metodi meshless dell ordine di convergenza e con ostacoli nel considerare discontinuit nelle propriet materiali e alcuni metodi non funzionano con distribuzioni arbitrarie di punti neces sitando di griglie regolari o richiedono specifici requisiti nella connettivit locale tra i punti per ottenere risultati precisi L aspetto dell onere computazionale nella definizione delle connettivit locali da taluni autori posto in particolare evidenza anche sotto lo specifico profi lo tassonomico della classificazione di un metodo come meshless se da un lato l idea centrale dell assenza di mesh uno degli elementi comunemente rite nuti focali Idelsohn 44 afferma che un metodo meshless un algoritmo in cui la definizione delle funzioni di forma dipende solamente dalla posizione dei nodi ciono nostante l esigenza di una rapida definizione della connettivit locale tra
50. tra il calore generato Gp all interno della sua regione tributaria e il calore uscente dalle sue superfici di contorno avendo calcolato cella per cella del com plesso locale tutti i termini di flusso Q attraverso le porzioni di superficie Af di contorno della regione tributaria del nodo polo P si ha gt Op Gp 4 3 ceI P dove il termine a primo membro espresso in funzione dei valori nodali di temperatura di tutti i nodi satellite S della costellazione Si ottiene una espressione nella forma Tp Ts E ili 4 4 Ts N c Scfr 30 85 4 L approccio meshless con il Metodo delle Celle I coefficienti della matrice a primo membro andranno a posizionarsi nella matrice fondamentale del sistema lineare risolvente alla riga corrisponden te all indice di nodo asociateo al nodo polo nella numerazione complessiva globale dei nodi le colonne in cui si posizioneranno i vari coefficienti sono quelle corrispondenti agli indici dei nodi della costellazione polo e satelliti nella medesima numerazione Analogamente accade per il valore a secon do membro che va a costituire posizionarsi nel termine noto del sistema risolvente alla riga associata al nodo polo Nota come precisato nell introduzione la definizione della regione tributaria attorno al nodo polo arbitraria Nel caso siano presenti sorgenti distribuite ne cessario assicurare una piena copertura del dominio con l insieme delle regioni tributari
51. tradizionali proiezioni radiografiche 10 1 3 La definizione di un modello Figura 1 1 Immagine CT pixel e voxel 1 3 1 1 La ricostruzione dell immagine La ricostruzione dell immagine quel processo con cui partendo dai dati grezzi per ciascuna delle giaciture di proiezione utilizzate si ricostruisce la distribuzio ne spaziale dell assorbimento subito dalla radiazione incidente nell area spazzata dai raggi di proiezione Esiste una precisa relazione matematica tra la distribuzione spaziale delle pro priet fisiche in un mezzo quali la densit l impedenza acustica la magnetizza zione e le misure di attenuazione di un segnale acustico o elettromagnetico Il processo di ricostruzione il percorso inverso attraverso il quale dalle misure di attenuazione si ottiene la distribuzione spaziale delle propriet fisiche del mezzo Nel caso in cui si usino i raggi X un fascio collimato passa attraverso l oggetto in esame ed esce dalla parte opposta attenuato a causa dell assorbimento apporta to dai tessuti presenti Rilevatori registrano i segnali in uscita e li convertono in proiezioni lineari se condo una direttrice ortogonale alla direzione del fascio della distribuzione del coefficiente di attenuazione lineare della sezione Noto l insieme di profili di at tenuazione lineari utilizzando opportuni algoritmi si riesce a risalire al profilo di attenuazione spaziale sulla sezione spazzata dal fascio de
52. u P in un pun to P all interno della cella c 60 3 2 Problemi di elasticita lineare tridimensionali Raccogliendo le componenti di spostamento di tutti i nodi di cella nel vettore Uh Uh Wh Ui Vi Wi Uj Ue Uj Wj Uq Vq Wq 3 40 le espressioni riportate nella 3 39 possono venire sinteticamente indicate come u P B n Mu dove essendo I una matrice identit 3x3 si inteso indicare I 0 0 I I 0 I 0 I I 0 0 I 1 21 0 SEGA 1 21 1 21 I 0 1 2I I 0 0 I 1 21 0 I 0 1 29 e B m I Be i 20 0 0 0 0 0 0 I 0 0 0 0 0 0 1 41 0 0 0 0 1 22 0 0 1 4 0 0 1 21 0 1 4 IE In I Ign Inc Ice IE Ip 3 2 2 Spostamenti deformazioni 3 2 2 1 Notazione vettoriale 1 o oc c6 c COCO 0 I 0 0 1 41 1 41 0 0 4I 1 41 0 O oo 0 1 41 1 41 0 0 1 41 3 41 O So el OO Oo 3 42 1 41 1 41 1 41 IC 3 43 Lo stato deformativo in una porzione di solido completamente descritto da 6 parametri tre di dilatazione secondo gli assi di riferimento e tre di scorrimento 61 3 Formulazione quadratica con il Metodo delle Celle y fra coppie di questi assi Questi sei parametri sono generalmente raccolti in una matrice detta tensore di deformazione Exa Ery Erz E Sie Eyy Eys 3 44 Ezg Ezy Ezz Questo tensore strettamente legato al campo degli spostamenti pi precisamente esso corrisponde alla parte simmetrica della matrice gradien
53. viscoelastici e dell eventuale insorgenza di de formazioni irreversibili Le condizioni di prova erano le seguenti temperatura 29 31 C umidit rela tiva 40 Il campione era tenuto costantemente avvolto in panni impregnati di soluzione salina I dati grezzi uscenti dalle griglie estensimetriche sono stati preliminarmente esa minati con il criterio di Chauvenet per l eliminazione di dati inattendibili che per sono risultati assenti La ripetibilit tra le misure risultata ottima cos co me il test di linearit delle deformazioni di ciascun estensimetro in funzione del carico applicato Sono state quindi calcolate le medie delle cinque ripetizioni di 140 5 2 Il test esplorativo misura per ciascun canale di ciascun estensimetro 5 2 2 Il dataset CT Il femore stato sottoposto a tomografia computerizzata scanner HiSpeed Ge neral Electric Co USA parametri di scansione 120 kVp 180 mA tipici di un normale esame clinico dopo essere stato scongelato ed immerso in acqua per prevenire artefatti di indurimento del fascio di radiazioni stato applicato un protocollo di scansione assiale con spessore del fascio 1 mm e distanza assiale tra scansioni successive pari a 2 mm nella regione prossimale dalla testa al piccolo trocantere e di 4 mm nella regione diafisaria e distale Il valore dimensionale del pixel risultato di 0 46 mm 5 2 3 Il modello FEM subject specific Il modello ad elementi finiti stato
54. 0 le immagini raffiguranti il campo di spostamenti ottenuto con una distribuzione casuale di punti la distribuzione dei punti stata ottenuta con variazioni di posizionamento casuali compatibili con la geometria applica te ad una griglia equispaziata di partenza Va osservato che variando la distanza media locale dei punti nelle costellazioni man mano create una condizione del genere non si presti a confronti numerici sistematici mancando un parametro di qualificazione attendibile della densit dei punti 4 4 Interpolazione locale con funzioni di base radiale Nel paragrafo precedente stata presentato un approccio meshless che muoven do da alcune considerazioni attorno all equazione di bilancio sulle regioni tri butarie associate ad ogni nodo costruisce un sistema algebrico risolvente senza richiedere la costruzione di una mesh nemmeno a livello locale fissata l attenzio ne su un nodo si individua un opportuno gruppo di nodi circostanti con questo 18manca qui una dimensione rappresentativa della dimensione di cella e la distanza me dia dei satelliti dal polo risultato essere un parametro con fluttuazioni eccessive per risultare attendibile in misura tanto pi marcata all aumentare della distorsione della distribuzione dei punti 111 4 L approccio meshless con il Metodo delle Celle Figura 4 19 Test della trave di Timoshenko distribuzione casuale di punti gene rata dalla griglia r
55. 0007 Tabella 4 1 Valore dell errore medio per i punti interni al dominio al variare del passo della griglia di punti Questi stessi valori riportati in scala logaritmica nel grafico di figura 4 8 evidenziano come l ordine di convergenza per le varie interpolazioni sia pari a 2 relativamente alla griglia regolare di punti 4 2 2 Sistemi elastici piani in regime di piccoli spostamenti Come caso del problema scalare sia dato un sistema piano di geometria e pro priret note su cui siano applicate delle forze e siano imposti dei vincoli agli spostamenti All interno e sul contorno di questa geometria venga disposto un conveniente numero di punti la cui posizione sia nota relativamente ad un rife rimento cartesiano globale xOy come illustrato in figura 4 3 Lo scopo che ci si 88 4 2 Meshless a celle locali con il CM log 6 T T T T T T 1 1 6 1 4 1 2 1 0 8 0 6 0 4 0 2 log lato Figura 4 8 Errore medio rispetto al passo di griglia per le implementazio ni con 4 5 e 8 nodi satellite prefigge quello di determinare le componenti del vettore spostamento in tutti i punti posizionati La procedura sar sostanzialmente la stessa illustrata nel precedente paragrafo 4 2 1 con la differenza che entro ogni cella andr determinata la forza agente su ogni porzione di superficie delimitante la regione tributaria del nodo polo in fun zione dei valori delle componenti di spostamento nodale Per chiarez
56. 03 13 2001 20 S A Marom and M J Linden Computer aided stress analysis of long bones utilizing computed tomography J Biomech 23 5 399 404 1990 21 C Zannoni R Mantovani and M Viceconti Material properties assignment to finite element models of bone structures a new method Med Eng Phys 20 10 735 40 1998 22 T V Nguyen J R Center and J A Eisman Femoral neck bone loss predicts fracture risk independent of baseline bmd J Bone Miner Res 20 7 1195 201 2005 23 E F Morgan H H Bayraktar O C Yeh S Majumdar A Burghardt and T M Keaveny Contribution of inter site variations in architecture to trabecular bone apparent yield strains J Biomech 37 9 1413 20 2004 24 Marco Viceconti and Fulvia Taddei Automatic generation of finite ele ment meshes from computed tomography data Critical reviews in biomedical engineering 31 1 2 27 72 2003 25 T McInerney and D Terzopoulos Topology adaptive deformable sur faces for medical image volume segmentation IEEE Trans Med Imaging 18 10 840 50 1999 180 26 L Mosekilde S M Bentzen G Ortoft and J Jorgensen The predictive value of quantitative computed tomography for vertebral body compressive strength and ash density Bone 10 6 465 70 1989 27 C Bergot V Bousson A Meunier M Laval Jeantet and J D Laredo Hip fracture risk and proximal femur geometry from dxa scans Osteoporosis International 13 7 542 550 July 2002
57. 2000 Figura 6 5 Diagramma della relazione tra HU e modulo di elasticit E espresso dalle relazioni di cui al paragrafo 6 1 3 7 00e 05 6 00e 05 5 00e 05 4 00e 05 3 00e 05 2 00e 05 1 00e 05 0 00 965 00 644 35 323 70 3 05 317 60 638 25 958 90 1279 55 1600 20 1920 85 Figura 6 6 Istogramma delle frequenze di HU nel dataset di uno dei femori Per la scelta del valore soglia HU ci si rifatti al campo di validit sperimenta le della definizione dell equazione di calibrazione densitometrica questa viene 158 6 2 Le criticita emerse infatti determinata in fase di acquisizione delle immagini CT inserendo in coda al ai segmenti ossei un supporto che alloggia tre soluzioni di idrossiapatite di titolo e dunque di densit nota La linearit rilevata tra unit HU e densit ha un intervallo di validit il cui limite inferiore pocr 50 mg cm questo valore corrispondente nei casi in oggetto ad un valore HU pari a 68 assunto come soglia dell equazione 6 5 Va rilevato come questa assunzione abbia a priori un margine di arbitrarie 400 I HE I I 300 i I l 200 i I Il 100 300 AN_e3 400 Figura 6 7 Valori di deformazione in alcuni estensimetri nel modello MCM in una delle prove al variare del valore di soglia HU dell equazione 6 5 t per contro a supporto di questa scelta due considerazioni contribuiscono ad avvalo
58. 5 575 6 5 4 log 5 7 54 1 0 8 0 6 0 4 0 2 0 log lato Figura 4 25 Test della trave di Timoshenko errore medio rispetto al passo di griglia in scala logaritmica con approccio CM meshless interpolazione RBF 130 Capitolo 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle Nel capitoli precedenti sono stati esplorati alcuni sviluppi basati sul Metodo delle Celle la formulazione quadratica per problemi scalari e vettoriali tridimensionali e alcuni approcci meshless Lo sviluppo e la verifica dell interpolazione quadratica con il CM ha confermato le potenzialit del Metodo in ambito tecnico applicativo gli acclarati vantaggi propri dalla formulazione discreta diretta alla base del CM manifesti ad esem pio nelle applicazioni di elasto dinamica e nei problemi di frattura sono disponi bili con una implementazione caratterizzata da precisioni di calcolo non inferiori a quelle raggiungibili dall affermato Metodo degli elementi Finiti Gli approcci meshless proposti aprono peraltro interessanti prospettive i vantag gi dell assenza di una fase di mesh o di remeshing appaiono accessibili senza le difficolt nell imposizione delle condizioni al contorno e senza i problemi di com plesist formale che affliggono molti degli approcci meshless di tipo differenziale Gli studi esplorativi circa le implementazioni meshless assumono particolare in teresse nell ambit
59. Celle semplice ed espressivo di distribuzione p dimensionale Essa l analogo di una funzione di campo usata nella formulazione differenziale la quale associa ad ogni punto non ad ogni p cella la densit q della grandezza Q Una distribuzione p dimensionale una funzione di dominio o funzione di insieme mentre una funzione di campo una funzione di punto 2 6 Impostazione del problema con il Metodo delle Celle Un problema fisico viene ben posto quando siano precisati e forma e dimensione del dominio e natura e posizione delle sorgenti e le caratteristiche materiali e le condizioni al contorno Risolvere il problema di campo significa allora determinare la distribuzione spa ziale e temporale della variabile principale di configurazione denominata in ge nere come potenziale del problema in esame Qualunque strada si segua per affrontare la soluzione del problema si perviene sempre alla formulazione di un equazione fondamentale algebrica o differenziale a seconda dell impostazio ne adottata Questa equazione nasce dalla composizione di due categorie di equazioni e equazioni di struttura e equazioni costitutive 2 6 1 Equazioni di struttura Le equazioni di struttura sono le equazioni che forniscono la struttura di una teoria e ne forniscono per cos dire le colonne portanti Esse sono dette anche equazioni topologiche e legano di fatto grandezze associate a celle di ordine n a celle di ordine n 1 il c
60. Il software di mappatura La struttura delle due versioni del codice BONEMAT la stessa e le differenze risiedono nel modo in cui partendo dai valori CT viene calcolato il valore medio HU assegnato a ciascun elemento della mesh La sua struttura di base del codice pu essere divisa in tre fasi fondamentali di seguito illustrate Valutazione del valore medio di HU e Per ciascun elemento della mesh il codice calcola un valore uniforme HU sulla base dei dati CT tramite la formula ut Sy HU amp y 2 dV n Jy HU r 5 t det J r s t dV 12 es dV Va dove V il volume dell elemento n x y z sono le coordinate nel riferimento CT r s t sono le coordinate locali nel sistema di riferimento dell elemento e J lo jacobiano della trasformazione tra i due sistemi di coordinate Gli integrali dell equazione 1 2 vengono valutati numericamente con la possibilit di poter disponibile presso il sito http www biomedtown org B3C_Building products bonemat decritta analiticamente nella documentazione riportata nelle pubblicazioni di riferimento 21 20 1 3 La definizione di un modello scegliere l ordine dell integrazione numerica Il valore HU x y z in un generico punto del dominio CT determinato con un interpolazione tri lineare tra i valori degli otto punti adiacenti Questo metodo computazionalmente pi oneroso ri spetto al precedente ma permette di stimare il valore di HU per ogni elemento con un
61. UNIVERSITA DEGLI STUDI DI TRIESTE Sede Amministrativa della Scuola di Dottorato di Ricerca Posto di Dottorato attivato grazie al contributo dell Istituto Ortopedico Rizzoli di Bologna XXII CICLO DELLA SCUOLA DI DOTTORATO DI RICERCA IN INGEGNERIA CIVILE E AMBIENTALE INDIRIZZO INGEGNERIA DELLE STRUTTURE SVILUPPO VALUTAZIONE E APPLICAZIONE DI METODI NUMERICI ALTERNATIVI AL METODO DEGLI ELEMENTI FINITI IN PROBLEMI DI BIOMECCANICA ORTOPEDICA Settore Scientifico Disciplinare ICA R 08 DOTTORANDO RESPONSABILE DOTTORATO DI RICERCA ING MARTINO PANI PROF IGINIO MARSON RELATORE PROF ENZO TONTI UNIVERSIT DEGLI STUDI DI TRIESTE ANNO ACCADEMICO 2008 2009 Indice Introduzione 1 La creazione di modelli numerici subject specific di segmenti ossei 1 1 1 2 1 3 Scopi e ragioni dell analisi issues G9 wae Lana Criteri di applicabilit ate ra e Lio i La definizione di un modello 1 3 1 La Tomografia Computerizzata CT 1 3 2 La costruzione della mesh dai dati CT 1 3 3 Definizione delle propriet meccaniche 1 3 4 Iproblemi dell approccio FEMmesh 2 Il Metodo delle Celle 2 1 2 2 2 9 2 4 2 5 2 6 Le variabili fisiche esi Pala a oS we hk Be 2 1 1 Costanti e variabili ooa TA eae oe OR 2 1 2 Classificazione in base alle entit geometriche di riferimento 2 1 3 La classificazione delle variabili in base al ruolo L orientazione degli elementi
62. a correlazione tra HU e densit del tessuto osseo individuata in fase di acquisizione della tomografia tramite lo disponibile su www biomedtown org 153 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle European Spine Phantom 82 risultata essere pocr 0000788343 HU 0 00422356 sl 6 1 essendo pgcr il valore di densit radiologica La relazione densit elasticit stata assunta dalla letteratura 83 E 6 950 p GPa 6 2 app dove Papp rappresenta la densit apparente del tessuto La relazione tra densit radiologica pgcr e densit minerale papp e tra questa e la densit apparente papp Sono state assunte dalla letteratura 87 PQCT 1 14 ash 0 09 Pash 9 6 Papp ee Componendo queste relazioni si ottiene E a h a HU b k con a 0 00069152894737 b 0 07524249122807 6 4 h 0 87719298245614 k 0 078947368421053 a 6850 B 1 49 Il modulo di Poisson stato assunto pari a 0 3 84 I modelli sono stati registrati spazialmente con il sistema di riferimento spe rimentale misurando e riportando la posizione di alcuni punti di riferimento estensimetri livello del cemento punto di applicazione del carico L acquisizio ne della posizione dei punti avvenuta tramite digitizer MicroScribe 3DX im mersion Corporation San Jose CA USA con risoluzione pari a 0 2 mm l errore di registrazione avvenuta con algoritmo c
63. a a 10 nodi ha dimensioni 30 x 30 Scrivendo il bilancio su tutte le N celle che descrivono la geometria in esame si ottengono altrettanti sistemi del tipo dell equazione 3 66 che opportunamente assemblati portano alla scrittura del sistema globale F KU 3 67 dove e Fil vettore delle forze applicate a tutti gli N nodi del sistema F 24s Figs ia dial 67 3 Formulazione quadratica con il Metodo delle Celle e K una matrice di dimensioni 3Nx3N ed detta matrice di rigidezza globale e U il vettore degli spostamenti nodali T U t1 U1 W1 tee UN UN WN 3 2 5 1 Orientazione delle porzioni di superficie duale Nel calcolo delle forze che agiscono sulla superficie di confine delle regioni tribu tarie di ciascun nodo fondamentale che cella per cella venga assunta un unica orientazione valida per entrambi i nodi separati da ogni porzione di superficie Questa orientazione d infatti il segno della forza ottenuta dal prodotto fra il tensore degli sforzi e il vettore area L orientazione in questione un orientazione esterna di un elemento la superfi cie che fa parte del complesso duale Come tale essa viene naturalmente indotta quando sia stata assunta un orientazione interna per gli elementi del complesso primale Provvedere ad una convenzionale orientazione interna degli elementi del complesso primale quindi operazione non formale bens sostanziale per la costruzione dell equazione di bi
64. a che separa due immagini adiacenti La seconda circostanza importante che i valori di HU del dataset variano da punto a punto Come illustrato nel dettaglio nel seguente paragrafo 5 2 5 le unit HU sono correlabili con la densit minerale e questa a sua volta correlabile con un valore di modulo di elasticit si ha quindi disponi bile una mappatura di valori di modulo di elasticit variabile da punto a punto sui punti di una griglia regolare Figura 5 1 Apparecchiatura per scansione CT la scansione dopo opportune rie laborazioni fornisce una pila di immagini registrate rispetto ad un sistema di riferimento disposto come indicato in figura 5 1 La selezione delle metodologia meshless applica ta Nel capitolo precedente sono stati proposti tre differenti approcci meshless Di questi solamente il primo quello basato sulla creazione di un complesso di celle 133 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle locale si verificato adatto ad un applicazione con le direzioni sopra specificate costruzione diretta del modello sulla griglia regolare di punti del dataset CT e l approccio con interpolazione polinomiale costringe nei problemi tridi mensionali all adozione di polinomi interpolanti che o producono una ma trice di interpolazione singolare quando i punti sono disposti su una griglia regolare o hanno un elevato grado parziale e totale con significativa
65. a dei solidi si distinguono dalle variabili locali per il fatto di non essere n densit n tassi 2 1 3 La classificazione delle variabili in base al ruolo Un altra importante classificazione avviene invece in relazione al ruolo che le va riabili hanno nel descrivere il campo fisico Sotto questo profilo si individuano variabili energetiche variabili di configurazione lavoro calore energia potenziale energia cinetica energia interna energia libera energia elettrom potenza densit di energia vettore di Poynting funzione di Lagrange funzione di Hamilton azione ece variabili geometriche e cinematiche coordinate generalizzate spostamento velocit tensore di deformazione velocit di deformazione vorticit vel angolare potenziale velocit acceleraz gravit frequenza vettore d onda fase differenza di fase potenziale gravitazionale temperatura affinit potenziale chimico potenziale elettrico vettore campo elettrico flusso magnetico induzione magnetica ecc parametri o costanti Figura 2 1 capacit termica costante dielettrica carica dell elettrone massa del protone conducibilit termica resistenza elettrica modulo di elasticit temperatura di fusione numero di Avogadro numero di Reynolds velocit del suono indice di rifrazione velocit della luce costante di Planck viscosit ecc variabili di sorgente variabili sta
66. a forma generale Tx di dy 2 dy 3 dy 4 dy 5 dy 6 Ex Oy dai d22 d23 do 4 d25 d 6 Ey Oz d31 d3 2 d3 3 d3 4 d3 5 d3 6 Ex 3 55 Tay dyi d4o d43 dad das da 6 Yay Tyz ds ds 9 d5 3 ds 4 d5 5 ds 6 Yyz Tra dea de 2 de 3 de de dee Yeu Per un materiale a comportamento elastico lineare e isotropo la matrice D assume la forma LU 1 1 0 0 0 1 1 0 0 0 E 1 1 0 0 0 De 1 v 1 2 0 0 0 5 0 0 Gon 0 0 0 0 0 0 0 0 0 0 3 2 Problemi di elasticita lineare tridimensionali dove E e v sono rispettivamente il modulo di elasticit del materiale ed il suo coefficiente di contrazione trasversale altrimenti detto modulo di Poisson 3 2 4 Legame forze sforzi In una regione caratterizzata da sforzo uniforme la forza di superficie che agisce su una superficie piana A descritta dal vettore area a 4 0 0 viene espressa dalla nota relazione di Cauchy che in forma discreta esprimibile come Ei C Oyy Oe X ay 3 57 T Oze Ozy Ozz 0z ossia in notazione sintetica T Ma 3 58 Adottando per o la notazione vettoriale si ha in modo del tutto equivalente Ox T az 0 0 a 0 a T VO ag 0O as a lt 0 My 3 59 To 0 0 a O ay ar Ay Taz Tea Nel caso trattato all interno di ciascuna cella il campo delle tensioni risulta essere funzione lineare delle coordinate In questo caso come noto si ha la notevole propriet LS ev ndA ta Ya za a 3 60 A prodotto fra la matrice calcolata nel
67. a risolvente Nulla assicura infatti il carattere mutuo del la relazione polo satellite un nodo che figura nella lista dei satelliti di un dato polo non necessariamente contempler questo polo nella lista dei suoi satelliti Non sussiste quindi necessariamente la simmetria di ruoli tra i nodi nella scrittura 2la geometria delle celle solitamente espressa da un parametro denominato aspect ratio de finito come il rapporto tra la lunghezza del lato pi lungo e quella del lato pi corto di una elemento celle molto distorte individuate da elevato valore diaspect ratio sono caratterizzate da una matrice di interpolazione mal condizionata questo fatto che verr approfondito nel seguen te paragrafo determina una perdita di accuratezza nell interpolazione della funzione incognita e dunque in ultima analisi un decadimento dell accuratezza della soluzione numerica a livello globale di intero sistema 3o le righe nel caso in cui la variabile di configurazione sia di natura vettoriale 81 4 L approccio meshless con il Metodo delle Celle dei relativi bilanci e questo si traduce evidentemente nella perdita di simmetria nella matrice che riassume i bilanci nodali Figura 4 2 Generazione delle mesh locali e definizione delle regioni tributarie i nodi primali sono indicati con il simbolo il contorno del dominio in linea trat teggiata in linea continua sono definite le celle locali e in linea spessa le celle duali
68. aedri della mesh locale di ciascuno dei due nodi su ognuna delle quattro porzioni di superficie si calcola la forza agente relativamente all interpolazione del campo di spostamenti dei tetraedri di riferimento appartenenti alle mesh locali dei due nodi in esame Questa procedura adottata e valutata nella sua applicazione nel prosieguo di questa tesi viene qui denominata Meshless Cell Method MCM Va rilevato come questa procedura non sia affatto semplice in termini algoritmici per una disposizione arbitraria di punti nello spazio Nel caso specifico dell ap plicazione per la fortunata circostanza dell allineamento dei punti rende il pro cesso di mediazione delle forze molto semplice Un altro elemento che risulta qui significativo sottolineare come questa proce dura conservi comunque la sua raison d tre nella filosofia meshless i tetraedri sono qui delle infrastrutture per l interpolazione locale del campo degli spo stamenti nella regione di spazio a loro associata in questa regione le forze di superficie sono calcolate con riferimento all interpolazione del corrispondente te traedro anche in punti ad esso esterni In questo senso si pu considerare l VMCM come una sorta di cross over tra la procedura a celle locali in cui si considera solamente la regione di spazio interna ad ogni cella e quella ad interpolazione polinomiale descritte nel capitolo 4 Replicando i test della trave di Timoshenko con le stesse modalit
69. alla scrit tura del sistema risolvente La prospettiva locale nella scrittura dell equazione di bilancio nodo a nodo suggerisce un approccio meshless reso immediatamente possibile dal Metodo delle Celle Si prenda infatti in esame un generico nodo P interno ad un domi nio bi o tridimensionale studiato in cui siano stati posizionati un certo numero di nodi ad esempio quelli generati da una tradizionale generazione di celle pri mali chiameremo convenzionalmente questo nodo preso in esame nodo polo Attorno a questo nodo possibile individuare un certo numero di nodi scelti tra i nodi prossimi in modo tale da circondare il nodo polo evitando cio che i punti individuati si addensino attorno ad una stessa direzione rispetto al polo Chiameremo convenzionalmente questi nodi nodi satellite Con questa costellazione dei nodi satellite opportunamente ordinati possi bile costruire un insieme di celle che e condividano tutte il nodo polo 80 4 2 Meshless a celle locali con il CM e siano tra loro adiacenti e coprano interamente l intorno del nodo polo senza lasciare porzioni di do minio scoperte L insieme di queste celle costituisce una mesh locale ovvero un complesso pri male locale associato al nodo polo Il numero delle celle locali ossia il numero di nodi satellite la loro forma e la loro dimensione sono elementi di principio arbitrari il numero di satelliti dipen de da un la
70. ana P Pastor J J Barrau S Laroze and J P Mo rucci Development of a three dimensional finite element model of a human tibia using experimental modal analysis J Biomech 24 6 371 83 1991 5 Jean Daniel Boissonnat and Bernhard Geiger Three dimensional recon struction of complex shapes based on the delaunay triangulation volume 1905 pages 964 975 SPIE 1993 6 William E Lorensen and Harvey E Cline Marching cubes A high re solution 3d surface construction algorithm SIGGRAPH Comput Graph 21 4 163 169 1987 7 D Taylor and J H Kuiper The prediction of stress fractures using a stressed volume concept J Orthop Res 19 5 919 26 2001 8 M J Ciarelli S A Goldstein J L Kuhn D D Cody and M B Brown Evaluation of orthogonal mechanical properties and density of human tra becular bone from the major metaphyseal regions with materials testing and computed tomography J Orthop Res 9 5 674 82 1991 9 R J McBroom W C Hayes W T Edwards R P Goldberg and A A White Prediction of vertebral body compressive fracture using quantitative computed tomography J Bone Joint Surg Am 67 8 1206 14 1985 10 T S Keller Predicting the compressive mechanical behavior of bone Journal of Biomechanics 27 10 1159 68 1994 11 D R Carter and W C Hayes The compressive behavior of bone as a two phase porous structure J Bone Joint Surg Am 59 7 954 62 1977 179 12 J C Ri
71. aso del gradiente di una funzione grandezza as sociata alla linea che congiunge due punti in cui sono noti dei valori scalari le equazioni circuitali associano invece ad una superfici una grandezza correlata a grandezze associate alle linee le equazioni di bilancio realizzano infine una rela zione fra grandezze associate a superfici con grandezze associate ai volumi Le caratteristiche di queste equazioni sono 39 2 Il Metodo delle Celle e legare tra di loro variabili dello stesso tipo variabili di configurazione con variabili di configurazione variabili di sorgente con variabili di sorgente e non contenere costanti materiali e essere valide a qualsiasi scala dimensionale La validit ad ogni scala dimensionale significa in particolare che le equazioni di struttura sono le stesse a prescindere dalla forma dalle dimensioni e dalla natura del mezzo materiale della porzione di campo sul quale sono scritte In particola re si ha che le equazioni di struttura valgono inalterate tanto per la trattazione differenziale quanto per quella discreta Equazioni di bilancio Le equazioni fondamentali di molte teorie fisiche sono basate su un equazione di bilancio Considerata una regione chiusa dello spazio V in cui sia definita una generica grandezza estensiva Q l equazione di bilancio esprime le variazioni del la grandezza Q che si realizzano all interno del volume V in un dato intervallo di tempo T Nella sua forma generale un
72. ata l applicazione della interpolazione con le Funzioni di Base Radiale RBF nell approccio meshless del Metodo delle Cel le applicato allo studio della deformazione di sistemi elastici piani in regime di piccoli spostamenti La procedura assolutamente simile a quella seguita utiliz zando l interpolazione polinoamiale Assunto un dominio di geometria e propriet modulo di elasticit E e coeffi ciente di Poisson v note si provveda a posizionarvi all interno un conveniente numero di punti di seguito denominati nodi La posizione di questi punti ver r riferita ad un sistema di riferimento cartesiano esterno XOY Preso in esame un nodo interno generico chiamato convenzionalmente nodo polo si selezionino Nc di nodi prossimi disposti attorno al nodo questi nodi selezionati prenderanno il nome di nodi satellite e assieme al nodo polo co stituiranno la costellazione locale C di punti con cui verranno approssimati local mente i campi scalari delle componenti del vettore spostamento attorno al nodo polo 123 4 L approccio meshless con il Metodo delle Celle 4 4 2 1 Interpolazione Le componenti del campo degli spostamenti costituiscono due campi scalari che verranno approssimati con funzioni di base radiale centrate nei nodi della co stellazione Adottando per semplicita di espressione formale e senza che questo comporti perdita di generalit un riferimento cartesiano locale xPpy ottenuto traslando
73. ational Journal for Numerical Methods in Engineering 48 11 1549 1570 2000 46 Hui Ping Wang and Dongdong Wang Efficient meshfree computation with fast treatment of essential boundary conditions for industrial applications Journal of Engineering Mechanics 135 10 1147 1154 2009 47 S R Idelsohn E Onate and F Del Pin The particle finite element method a powerful tool to solve incompressible flows with free surfaces and breaking waves International Journal for Numerical Methods in Engineering 61 7 964 989 2004 48 Sergio R Idelsohn and Eugenio Onate To mesh or not to mesh that is the question Computer Methods in Applied Mechanics and Engineering 195 37 40 4681 4696 2006 John H Argyris Memorial Issue Part I 49 Vinh Phu Nguyen Timon Rabczuk Stephane Bordas and Marc Du flot Meshless methods A review and computer implementation aspects Mathematics and Computers in Simulation 79 3 763 813 2008 50 Yi Je Lim and Suvranu De Real time simulation of nonlinear tissue respon se in virtual surgery using the point collocation based method of finite sphe res Computer Methods in Applied Mechanics and Engineering 196 31 32 3011 3024 2007 Computational Bioengineering 51 Suvranu De Jung Kim Yi Je Lim and Mandayam A Srinivasan The point collocation based method of finite spheres pcmfs for real time surgery si mulation Computers amp Structures 83 17 18 1515 1525 2005 Advances in
74. ativa dell acquisizione CT pixel pari a 0 57um lad dove il modello FEM costruito con elementi di dimensione media dell ordine del millimetro Lo studio ha individuato dunque un processo di modellazione numerica di segmenti ossei caratterizzato da un elevato grado di automaticit il modello nu merico viene costruito direttamente sul dataset CT in maniera indipendente dallo specifico segmento osseo rappresentato la procedura non richiede alcuna mani polazione e si dimostrato poter ottenere accuratezze nella replicazione delle misure sperimentali in vitro assolutamente concordi con quelle della procedura FEM di riferimento Sebbene la dimensione computazione dei modelli sia signi ficativamente maggiore di quella dei corrispondenti modelli FEM in virta del l utilizzo della risoluzione nativa dei dataset CT e sebbene quindi il tempo di calcolo di ogni singolo modello sia dell ordine delle ore va rilevato come l eleva ta automaticit del processo rappresenti comunque un progresso relativamente all applicazione su vasta scala l impostazione di un calcolo richiede infatti po che elementari operazioni selezione della regione di interesse e imposizione di carichi e vincoli che necessitano di pochi minuti di lavoro di un operatore essen do il resto del processo costruzione assemblaggio e soluzione del sistema lineare risolvente completamente autonomo per contro un analisi FEM computazinal mente molto pi leggera e di calcolo molt
75. avoro invece piuttosto una riflessione concentrata sulle alternative al FEM che presentino carattere di innovazione Sotto questo aspetto accettando tem poraneamente la denominazione di Metodo delle Celle in attesa che la comunit 10come riportato nella presentazione del lavoro la valutazione che ha portato a identificare nel FEM il metodo numerico di elezione per lo studio di problemi di biomeccanica ortopedica ha nel passato affrontato una valutazione critica di tutte le alternative disponibili tra quelle consolidate per la categoria di problemi affini 95 4 L approccio meshless con il Metodo delle Celle scientifica affronti e risolva la questione gli approcci proposti rappresentano si curo motivo di interesse tanto da giustificare la loro presenza in questo elaborato esplorativo Nel seguito verranno quindi presentate due formulazioni che implementano la procedura descritta La prima delle due formulazione adotta come funzione in terpolante i polinomi essa presentata in 60 relativamente alla soluzione dell e quazione di Laplace in domini bidimensionalli viene qui riportata sviluppata e applicata alla soluzione del problema della deformazione elastica di sistemi pia ni in regime di piccoli spostamenti di seguito denominato per mere ragioni di sintesi come problema dell elasticita piana Il secondo approccio muove dall intenzione di superare alcuni limiti dovuti al l interpolazione polinomiale e
76. avvenuta replicando un test meccanico condotto in vitro su un femore umano con un modello numerico costruito direttamente da un dataset CT il modello re gistrato spazialmente rispetto al riferimento sperimentale replicava le condizioni di carico e di vincolo del test meccanico la valutazione avvenuta per compa razione dei valori di deformazione registrati sperimentalmente in punti anato micamente significativi con i valori predetti dal modello numerico nei medesimi punti Un ulteriore elemento di valutazione stato il confronto con le predizio ni di un analogo modello FEM replicante il medesimo esperimento ed eseguito secondo una procedura standard precedentemente validata rispetto alle misure sperimentali In questo primo test emersa una buona corrispondenza delle predizioni del mo dello MCM con i valori sperimentali gli elementi usati come valutazione coef ficiente angolare e intercetta della retta di regressione dati numerici rispetto ai valori sperimentali errore medio e massimo delle predizioni numeriche rispetto ai valori misurati sono risultati confrontabili con quelli caratterizzanti il modello FEM validato Questa prima valutazione risultata incoraggiante l ipotesi che l MCM possa co stituire un processo di modellazione automatico e accurato potenzialmente utile per la predizione dello stato di sollecitazione in segmenti ossei in vivo Bench il confronto abbia riguardato 13 differenti siti di misura delle deformaz
77. bbe grandi sforzi e risulterebbe eccessivamente onerosa in termini di tempo Un generatore automatico di mesh AMG automatic mesh generator pu essere definito 1 come un algoritmo automatico che pu creare una mesh valida per una geometria di complessit arbitraria con il solo controllo di parametri di pro gramma e senza l intervento manuale dell utente Ogni AMG pu essere classificato in relazione al tipo di connettivit della mesh strutturata o non strutturata 2 e in relazione al tipo di elemento utilizzato per la discretizzazione nella maggior parte dei casi tetraedri o esaedri Pur essen do un tema ancora aperto sono oggi disponibili per mote applicazioni algoritmi assai robusti per la generazione di mesh e quasi tutti i codici FEM commercia li permettono di generare mesh di tetraedri In particolari contesti applicativi quale per l appunto quello biomedico esigenze e problemi specifici limitano di fatto l applicabilit degli strumenti disponibili per la discretizzazione automatica Nella specificit della rappresentazione di un segmento osseo i maggiori ostacoli sono costituiti dalla complessit geometrica e dalla scarsa definizione con cui descritta la superficie di contorno della struttura ossea La CT oggi la fonte di informazioni pi accurata per la rappresentazione di una Hounsfield fu l ingegnere britannico che mise a punto la tecnologia della TAC operazione che gli valse il premio nobel per la med
78. calmente i valori delle variabili di configurazio ne ai nodi della costellazione sono e stabilit del problema di interpolazione rispetto alla disposizione spaziale dei punti della costellazione per un opportuna scelta delle funzioni di base radiale e indipendenza dal numero di nodi che compone la costellazione locale la funzione interpolante somma delle singole funzioni di base radiale cen trate nei vari punti della costellazione questa natura additiva della inter polazione rende assolutamente libero e variabile di volta in volta il numero dei nodi scelti per comporre la costellazione locale e indipendenza dalle dimensioni del problema 9 cfr 67 20cfr 68 115 4 L approccio meshless con il Metodo delle Celle L idea di usare le RBF per l interpolazione locale viene da quanto proposto in let teratura per metodologie meshless differenziali 69 70 71 molta letteratura sull argomento oltre che in ambito speculativo circa le loro propriet matema tiche 72 73 viene da ambiti applicativi quale la manipolazione di immagini 74 e lo sviluppo di reti neurali 75 dove alcuni aspetti funzionali vengono af frontati con la concettualizzazione di problema di curve fitting Nei paragrafi seguenti verranno illustrate le applicazioni dell interpolazione con RBF all implementazione meshless del metodo delle celle nella sostanza la proce dura del tutto identica a quella presentata per l interp
79. ce S C Cowin and J A Bowman On the dependence of the ela sticity and strength of cancellous bone on apparent density J Biomech 21 2 155 68 1988 13 Dianna D Cody Fu J Hou George W Divine and David P Fyhrie Femo ral structure and stiffness in patients with femoral neck fracture Journal of Orthopaedic Research 18 3 443 448 2000 14 K J Bozic J H Keyak H B Skinner H U Bueff and D S Bradford Three dimensional finite element modeling of a cervical vertebra an investigation of burst fracture mechanism J Spinal Disord 7 2 102 10 1994 15 J C Lotz E J Cheal and W C Hayes Fracture prediction for the proxi mal femur using finite element models Part i linear analysis Journal Of Biomechanical Engineering 113 4 353 360 1991 16 M Dalstra R Huiskes and L van Erning Development and validation of a three dimensional finite element model of the pelvic bone Journal Of Biomechanical Engineering 117 3 272 278 1995 17 A A Edidin Modeling the distribution of bone moduli in bone implant systems Cornell University Ithaca NY 1991 18 B Merz P Niederer R Muller and P Ruegsegger Automated finite element analysis of excised human femora based on precision qct J Biomech Eng 118 3 387 90 1996 19 P M Cattaneo M Dalstra and L H Frich A three dimensional finite ele ment model from computed tomography data a semi automated method Proc Inst Mech Eng H 215 2 2
80. cettuale tra la filosofia dell approccio discreto e la natura dell infor mazione usata per la costruzione dei modelli numerici in ultima analisi il CM aveva manifestato di poter rendere accessibili alcuni strumenti quali la formula zione quadratica con ordine di convergenza nativamente pari a 4 e un approccio meshless di semplice formulazione tali da costituire a priori motivo di interesse e ragione di studio in un ottica non gi solamente accademica quanto piuttosto pratica e applicativa Lo studio ha dapprima sviluppato l estensione della formulazione quadratica ai problemi scalari e vettoriali tridimensionali la sua applicazione alla soluzione numerica dell equazione di Laplace e alla soluzione di problemi di elasticit li neare ha confermato quanto gi registrato nel caso bidimensionale relativamente ai problemi scalari un opportuna costruzione del complesso duale consente al CM di di ottenere nativamente ordini di convergenza superiori a quanto teori camente previsto per il FEM nella sua formulazione ordinaria ciononostante in entrambi i casi le soluzioni sono risultate del tutto confrontabili con quelle ottenu te da implementazioni FEM di codici commerciali avanzati rendendo i risultati significativi dal punto di vista teorico ma non utili sotto il profilo applicativo in mancanza una maggiore precisione numerica ottenibile a parit di dimensio ne computazionale del problema prevalgono i vantaggi costituiti da un metodo di con
81. che completano la defi nizione del sistema fisico in esame Questa la manifestazione di una stretta corrispondenza che sussiste fra la classificazione fisica delle variabili e la classificazione geometrica La sistematica associazione fra variabili fisiche ed elementi spaziali orientati dei due sistemi di celle un punto fondamentale per la formulazione finita diretta delle leggi fisiche dei campi Resta ora evidente una sostanziale differenza fra i complessi di celle e le me sh tradizionali e costruire un complesso di celle significa certo dividere un dominio in tante porzioni discrete ma significa anche costruire il sistema di elementi geome trici duali a quelli appena creati Questa operazione che funzionale alla formulazione discreta usta dal Metodo delle Celle non ha un omologo nella pratica usuale dei metodi numerici citati e le variabili globali sono associate alle celle 0 celle vertici 1 celle lati 2 celle facce 3 celle volumetti Data una grandezza Q associata ad elementi spaziali di dimensione p ne viene che ad ogni cella p dimensionale cx del complesso di celle risulta associato un valore Q della grandezza Q In questo modo la distribuzione della grandezza Q descritta dalla N pla Q1 Qo Qy Nella topologia algebrica una tale distribuzione porta il nome quanto mai infelice di co catena Noi la denoteremo con un termine pi 38 2 6 Impostazione del problema con il Metodo delle
82. cipali di deformazione predetti dal modello MCM manifestano un e levata correlazione con i valori sperimentali R 0 85 La pendenza della retta di regressione valori numerici vs sperimentali risultata maggiore dell unit ma il valore dell intercetta risultato significativamente prossimo allo 0 come riportato e illustrato nel grafico di figura 5 11 dove per comparazione sono ri portati anche i valori ottenuti con il modello FEM e la relativa retta di regressione Gli indicatori di accuratezza globale sono risultati essere soddisfacenti essendo il valore di errore massimo pari al 40 e l errore quadratico medio apri al 16 co me riportati in tabella 5 1 Nella tabella 5 2 sono riportate le componenti dello spostamento estratte dai due modelli MCM e FEM nei due punti di controllo i valori suggeriscono un comportamento analogo dei due modelli avendo gli spo stamenti predetti la stessa direzione e ampiezze paragonabili A questo livello il modello MCM appare sistematicamente pi rigido del corrispondente modello FEM 146 5 3 Risultati 72000 y 1 10x 21 44 1500 FEM R 0 91 M y 1 21x 35 58 C R 0 85 e 2000 1500 1000 500 500 1000 1500 2000 2000 Figura 5 11 Retta di regressione dei valori di deformazione predetti ue asse y rispetto ai valori sperimentali ue asse x per MCM e FEM strain correlation MCM FEM R 0 85 0 91 pendenza 1 21 1 10 intercetta 36 21
83. co Schileo Enrico Dall Ara Fulvia Taddei Andrea Malandrino Tom Schotkamp Massimiliano Baleani and Marco Viceconti An accurate esti mation of bone density improves the accuracy of subject specific finite element models Journal of Biomechanics 41 11 2483 2491 2008 P J Besl and N D McKay A method for registration of 3 d shapes IEEE Transactions on Pattern Analysis and Machine Intelligence 14 239 256 1992 185
84. com plessit algebrica della formulazione e per distribuzioni regolari dei punti asimmetria nella distribuzione dei satelliti rispetto al polo e l approccio con interpolazione RBF certamente pi agile in termini di sele zione della costellazione risultato avere una complessit formale ecces siva per gli strumenti di manipolazione simbolica disponibili relativamen te al calcolo dell espressione analitica della equazione di bilancio al nodo in altri termini l espressione algebrica dell integrazione delle tensioni sul la superficie di una regione tributaria anche di forma elementare risulta ta eccessivamente complessa per risorse algoritmiche e computazionali dei manipolatori simbolici disponibili L approccio meshless a celle locali dunque il candidato per l applicazione di test Nel caso tridmensionale la mesh locale costruita attorno ad ogni nodo composta da tetraedri aventi vertice nel nodo polo e gli altri tre vertici nei nodi primi vicini La disposizione regolare dei punti sulla griglia strutturata del da taset CT fa assumere a questa mesh una conformazione di semplice definizione come illustrato in figura 5 2 L applicazione di questa metodologia ha per individuato problemi numerici legati alla non simmetria della matrice di rigidezza del sistema risolvente questa asimmetria origina dalla disomogeneit delle propriet dei materiali associate al le celle locali di ciascun nodo alle celle della mesh local
85. dalla adozione delle Funzioni di Base Radiale RBF per l interpolazione locale delle variabili di configurazione Le RBF introdotte per risolvere problemi di interpolazione multidimensionali costruiscono la funzione approssimante come combinazione lineare di un certo numero di funzioni che hanno come argomento reale positivo Si immagini che della funzione incognita f R R 4 59 siano noti solo N valori in altrettanti punti del suo dominio Si desidera individuare una funzione y definita in R a valori reali che approssimi la funzione incognita f assumendo i valori noti nei punti dati y R R ylz dp Peet 4 61 113 4 L approccio meshless con il Metodo delle Celle L idea delle RBF quella di costruire la funzione approssimante y come som ma di N valori assunti da una funzione reale di argomento la norma x xill paS gt wid Ile zll 4 62 dove e w R sono detti pesi e x punti in cui sono noti gli N valori d sono detti centri e d r una funzione a valori reali definita in R Circa la funzione r va precisato che le funzioni solitamente utilizzate sono e la funzione spline r r con n dispari la funzione thin plate spline r r logr con n pari e la funzione gaussiana r2 r e 202 a gt Q e la funzione multiquadratica inversa 1 1 2 b r c gt Q la funzione multiquadratica r 72 02 o D La funzione gaussiana e la fu
86. dell estensi metro saranno quindi ottenuti mediando i valori principali calcolati ai centroidi dei voxel dell insieme Yp Con la procedura descritta si in grado di ottiene un informazione che medi i va lori alla scala dell acquisizione considerando i soli voxel interni al tessuto osseo Assumendo R apri a 2 mm si riesce infatti a tenere in considerazione da un lato l estensione della superficie di acquisizione sperimentale della deformazione l e stensimetro ha forma circolare con diametro di circa 3 mm dall altro l incertezza del suo posizionamento nel riferimento del datset CT Il grafico di figura 6 8 illustra come varino le predizioni dei valori di deforma zione ai vari estensimetri variando il valore di R nel grafico vengono riportate 5 colonne per ogni misura la prima corrisponde al valore sperimentale la seconda alla predizione FEM la terza alla predizione MCM puntuale le ultime due alle predizioni MCM mediate con valore di R pari a 2 e a 3 mm si nota come la me diazione con la procedura descritta elimini la presenza di valori non significativi corrispondenti a informazioni puntuali lette nella matrice acquosa al contem po le predizioni con la procedura mediata risultano in prima approssimazione stabili alla variazione di R 6 3 Validazione con le misure sperimentali Con la procedura appena descritta i valori principali di deformazione predetti dal modello MCM per varie condizioni di carico nei vari femori hanno u
87. della trave di Timoshenko errore medio rispetto al passo di griglia in scala logaritmica con approccio CM meshless a celle locali Il concetto di bilancio di nodo condotto sulla regione tributaria costituita dalla corrispondente cella duale suggerisce infatti la definizione di una mesh locale che circondi ogni singolo nodo allo scopo di definirne l associata cella duale sua regione tributaria La definizione delle grandezze coinvolte nel bilancio di nodo avviene cella per cella con la consueta interpolazione delle grandezze nodali di cella Sotto questo profilo le celle primali del complesso primale locale definito attorno al nodo polo costituiscono quindi una infrastruttura per collocare e de finire le grandezze coinvolte nel bilancio locale del nodo polo In particolare le 94 4 3 Interpolazione locale con funzioni polinomiali singole celle sono il supporto per la interpolazione locale dei valori nodali delle ricercate variabili di configurazione e delle grandezze da esse derivate Questa considerazione unitamente all ispirazione offerta da alcune metodologie meshless in ambito prettamente differenziale suggeriscono la seguente procedu ra e fissata l attenzione su un nodo il nodo polo individuare un certo numero di nodi prossimi circostanti i nodi satellite definendo una costellazione di punti e interpolare i valori incogniti delle variabili di configurazione nei vari punti della costellazione con una
88. di orientazione interna gli enti duali di orientazione esterna Una grandezza associata agli intervalli che cambi segno per inversione tempo rale necessita evidentemente di un orientazione interna dell intervallo di tempo essendo legata al verso di percorrenza Peraltro non si vede il significato fisico del cambiamento del verso per un intervallo dotato di orientazione esterna La velocit che cambia segno invertendo la direzione di percorrenze dell asse dei tempi viene indubbiamente associata a intervalli primali mentre l impulso in variante rispetto al verso di percorrenza legato ad intervalli duali 35 2 Il Metodo delle Celle In conclusione si pu rilevare come le variabili di configurazione risultino si stematicamente associate elementi dotati di orientazione interna per contro le variabili di sorgente risultano associate ad elementi dotati di orientazione ester na 2 5 complessi di celle Nella trattazione differenziale si utilizzano le funzioni di punto e risultano per tanto fondamentali i sistemi di coordinate utili a collocare i punti nello spazio e quindi a collocare le grandezze fisiche La trattazione discreta considera invece grandezze globali che sono come det to riferite a punti linee superfici e volumi Il riferimento spaziale costituito dai sistemi di coordinate risulta pertanto uno strumento insufficiente necessaria invece una struttura organizzata di elementi geometrici volumi facce spi
89. dio rispetto al lato della mesh per CM e FEM y 1 88x 1 38 R 0 999 log 5 bg 2 6 4 a FEM 3 4 a x CM 3 2 al T T T T 17 0 6 log I 0 5 0 3 0 2 Figura 3 11 Sfera cava pressurizzata grafico errore lato l in scala logaritmica 71 3 Formulazione quadratica con il Metodo delle Celle 72 Capitolo 4 L approccio meshless con il Metodo delle Celle In questo capitolo vengono formulate alcune ipotesi di approccio meshless del CM sia epr i problemi scalari che per problemi di elasticita lineare in domini bi dimensionali La prima parte del capitolo dedicata ad una breve introduzione sui metodi me shless in ambito differenziale conclusa da una digressione sulle applicazioni nel lo specifico ambito della biomeccanica I caratteri generali di un approccio me shless reso accessibile in modo diretto dalla struttura concettuale del CM ven gono poi declinati nello specifico di tre approcci proposti il primo basato sulla creazione di un complesso locale di celle primali gli altri due basati su un in terpolazione della funzione incognita definita localmente su una costellazione di punti a connettivit non definita 4 1 La filosofia dei metodi meshless Il Metodo degli Elementi Finiti FEM si ormai affermato come il principale e pi diffuso strumento di analisi numerica per soluzione di una vastit di pro blemi scientifici e tecnico ingegneristici un tempo inaffrontab
90. discretizzazione del segmento osseo necessario attri buire ad ogni elemento della mesh le relative propriet meccaniche Nel caso di una modellizzazione generica di un osso le propriet meccaniche dei diffe renti tessuti ossei sono generalmente assunti pari ai valori medi riportati nella letteratura delle misure sperimentali 4 7 Quando invece si cura un model lo subject specific le propriet meccaniche devono essere derivate dai dati CT stato dimostrato come la distribuzione delle tensioni in un osso sia fortemente legata alla distribuzione delle propriet meccaniche nel tessuto osseo Pertanto di fondamentale importanza individuare un metodo valido per individuare dai dati CT la distribuzione delle propriet meccaniche in un tessuto osseo e poterla mappare in un modello ad elementi finiti subject specific In una CT a raggi X le immagini rappresentano la distribuzione del coefficien te di attenuazione lineare dei tessuti i valori riportati sono in genere espressi in Hounsfield Units HU unit di misura ormai riconosciuta come standard Que sto sistema di unit presenta i coefficienti di attenuazione nella forma normaliz zata definita dall espressione 1 1 nella quale all acqua viene assegnato il valore 0 e all aria il valore 1000 stato dimostrato che esiste una correlazione quasi lineare fra le unit HU CT numbers e i valori di densit apparente di un tessuto biologico e che questo le game valido per u
91. do la forma debole In particolare tra le formulazioni meshless che usano il metodo di Galerkin la discretizzazione della forma differenziale pu avvenire o usando dei punti di in tegrazione in cui le funzioni sono valutate cosiddetti stress points oppure ricor rendo alla cosiddetta background mesh o cell structure un set di celle di integrazione ove viene utilizzato il metodo della quadratura di Gauss 75 4 L approccio meshless con il Metodo delle Celle Sebbene l assenza della mesh rappresenti oltre che ragione di novit anche motivo di attrattiva verso questi metodi va ciononostante rilevato come essi non si siano ancora affermati nella pratica corrente pochi sono i codici commerciali che implementano algoritmi meshless e molti metodi sono ancora in fase di svi luppo con applicazioni sperimentale in casi tecnologici la cui coplessit giustifichi lo sforzo di una consistente indagine teorica e speculativa La maggior parte dei metodi meshless sinora sviluppati accusa infatti alcuni pun ti critici e difficolt nell imporre le condizioni al contorno di tipo essenziale il carat tere non interpolatorio delle funzioni di forma rendono impossibile l impo sizione diretta di questo tipo di condizioni al contorno per cui si deve ricor rere a tecniche dedicate le tecniche pi diffuse sono la tecnica dei moltipli catori di Lagrange e i penalty methods ma va sottolineato come il tema sia ancora oggetto di indagine e di
92. domi nio studiato Per questa ragione quando il nodo polo giaccia sul bordo la sua re gione tributaria non lo circondera completamente ma solamente per quella parte che risulta interna al dominio figura 4 14 4 3 2 5 Testsulla trave di Timoshenko L indagine sulla precisione dell implementazione stata condotta relativamente al problema della trave di Timoshenko per il quale si dispone di una soluzione analitica di confronto 109 4 L approccio meshless con il Metodo delle Celle Due differenti test sono stati condotti relativamente allo stesso problema utiliz zando una costellazione locale composta da 4 e da 8 satelliti ossia una costella zione composta da 5 o da 9 nodi Geometria propriet materiali e carico usati nel test numerico sono riportati in tabella 4 5 per riprodurre con esattezza le condizioni di validit della soluzione teorica su tutti i punti del contorno sono stati imposti i valori di spostamento previsti dalle equazioni 4 13 nella posizione corrispondente Come metrica di precisione stato considerato il valore quadratico medio del modulo della dif ferenza vettoriale tra spostamento numerico e spostamento teorico in ciascuno dei punti non vincolati uisnum Wit Vinum vie 4 58 N avendo indicato con e N il numero dei nodi interni dunque non vincolati e per cui ha senso il confronto con la soluzione teorica e con u e con v le componenti del vettor
93. e formazione in un punto avviene derivando opportunamente la funzione interpo lante gli spostamenti nodali del voxel che contiene quel punto I valori principali 155 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle di deformazione usati per il confronto con i dati sperimentali sono quindi otte nuti calcolando i valori principali nei punti sede degli estensimetri e escludendo dei tre valori quello pi prossimo allo zero 6 2 Le criticit emerse Il confronto dei valori di ottenuti con i modelli MCM hanno evidenziato due criticit importanti della peocedura di modellazione descritta la prima legata all attribuzione delle propriet meccaniche desunte dai valori HU del datset la seconda di carattere pi specificamente interpretativo legata all estrazione delle informazioni per la valutazione di rispondenza con i dati sperimentali 6 2 1 Attribuzione delle propriet meccaniche nell MCM La comparazione macroscopica del campo di spostamenti nei modelli omologhi MCM e FEM evidenziava una sistematica e consistente maggiore rigidezza del modello MCM il quale registrava spostamenti massimi sistematicamente inferio ri al FEM Questo fatto era confermato anche dall analisi puntuale dei valori di deformazione sia in relazione ai dati sperimentali che in relazione alle predizioni FEM La causa di questo fenomeno stata individuata nella combinazione dell equa zione di correlazione tra HU e modu
94. e 0 cella duale in ogni 2 o 3 cella primale quanto il modo di unire tali punti collegamento diretto oppure appoggiato al baricentro 37 2 Il Metodo delle Celle dei lati o delle facce delle celle primali sono assolutamente arbitrari Le pi sem plici figure geometriche che compongono un complesso di celle nel piano sono i triangoli un complesso del genere viene pertanto detto complesso simpliciale il complesso simpliciale nello spazio invece costituito da tetraedri Per entram bi questi complessi la costruzione delle celle duali viene comunemente effettuata adottando come punti duali i baricentri oppure i circocentri sferocentri delle cel le primali Rispetto ai termini dell orientazione viene chiamato primale il complesso di cel le cui si attribuisce orientazione interna conseguentemente il complesso duale quello che riceve dal complesso primale un orientazione sterna Un altro aspetto della dualit dei sistemi di celle riguarda l attribuzione del le grandezze fisiche agli elementi geometrici se associamo alle celle del sistema primale le variabili fisiche di configurazione e attribuiamo dunque al complesso primale orientazione interna si produce una naturale associazione delle variabili di sorgente che come detto necessitano di orientazione esterna ai corrisponden ti elementi duali La costruzione del sistema duale mette cio a disposizione le strutture necessarie ad ospitare le variabili di sorgente
95. e di tetraedri costruita attorno ad un nodo va infatti associato un valore di modulo di elasticit che sia in qualche modo legato ai valori dei nodi su cui la mesh locale si appoggia qualunque strategia si adotti per individuare questo valore media dei valori ai nodi della costellazione media pesata valore associato al nodo polo si ha comunque che in ogni mesh locale sono presenti valori di E differenti da quel li associati alle celle delle mesh locali adiacenti anche se la griglia regolare e anche se dunque la relazione polo satellite risulta simmetrica le forze di su perficie calcolate sull elemento di separazione tra due nodi affiancati diversa quando venga calcolata relativamente alla mesh dell uno o dell altro nodo fatto sono stati utilizzati i programmi open source di Computer Algebra System CAS Ma xima http maxima sourceforge net Mathomatic http www mathomatic org math e il pacchetto Symbolic della suite open source Octave Forge http octave sourceforge net 134 5 1 La selezione delle metodologia meshless applicata Figura 5 2 Mesh locale costituita da tetraedri a 4 nodi appoggiati al nodo polo e ai nodi primi vicini questo che produce un asimmetria nella matrice di rigidezza globale Si verificato che realizzando una mediazione delle forze agenti sulla superficie di interfaccia tra due nodi vicini considerata appartenente alle regioni tributarie dei due nodi
96. e distribuite forze di volume infatti necessario che l insieme delle regioni tributarie definite localmente copra l intero dominio in modo che nessuna sua porzione e dunque nessun termine sorgente venga tralasciata Questo aspetto non banale quando le regioni tribu tarie di nodo siano definite dagli elementi locali di un complesso primale locale L interpolazione polinomiale locale invece usa una regione tributaria svincolata dalla disposizione della costellazione di punti individuata in particolare essendo completamente libera l estensione della regione tributaria la copertura dell inte ro dominio controllabile A fronte di questi vantaggi l interpolazione locale dei valori nodali con una funzione polinomiale ha per due fondamentali inconvenienti e il numero di nodi che compongono la costellazione fisso dipendendo dal polinomio adottato e non pu venire variato che per intervalli fissa ti si perde in questo modo la possibilit di migliorare l approssimazione arricchendo la costellazione dei nodi in dipendenza di fattori locali e la disposizione geometrica dei punti scelti non sempre compatibile con una accettabile accuratezza della interpolazione la matrice di interpolazio ne o di collocazione pu risultare singolare per particolari disposizioni dei nodi della costellazione rendendo laboriosa la scelta dei punti adatti Rispetto a queste limitazioni una possibile soluzione proviene
97. e ed intercetta 145 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle della retta di regressione individuata in condizioni ideali di accordo tra modello numerico e dati sperimentali vi dovrebbe essere perfetta linearit nella relazione R 1 con pendenza pari ad uno equivalente qui ad una corretta stima della rigidezza ed intercetta pari a zero assenza di errore sistematico di sovra o sottostima dei dati sperimentali e vengono calcolati alcuni parametri derivati dalle stime puntuali di errore nella predizione l errore quadratico medio RMSE il suo valore percen tuale rispetto al massimo valore di misura sperimentale l errore massimo e analogamente il suo valore percentuale rispetto al massimo valore misurato Come ulteriore elemento di valutazione volto ad delineare le differenze di comportamento tra MCM e FEM sono stati confrontati i valori dello spostamento calcolati dai due modelli in due punti anatomicamente ben identificabili il centro della testa del femore e il punto posteriore del piccolo trocantere 5 3 Risultati Verificato che non sussistesse dipendenza tra la condizione di carico e le diffe renze tra valori di deformazione sperimentali e quelli predetti dal modello MCM ANOVA fattoriale con P value gt 0 9 i valori di deformazione sono strati consi derati in maniera aggregata per tutte le condizioni di carico e analizzati nel loro complesso I valori prin
98. e equazioni 4 13 nella posizione corrispondente Come metrica di precisione stato considerato il valore quadratico medio del modulo della dif ferenza vettoriale tra spostamento numerico e spostamento teorico in ciascuno dei punti non vincolati 6 _ Cani fa ieee uri Vit 4 103 avendo indicato con e N il numero dei nodi interni dunque non vincolati e per cui ha senso il confronto con la soluzione teorica e con u e con v le componenti del vettore spostamento secondo la direzione x e y rispettivamente del sistema di riferimento globale e con il pedice num un valore di soluzione numerica e con il pedice un valore di soluzione teorica Il calcolo stato condotto considerando costellazioni locali composte da 8 nodi e utilizzando una distribuzione di punti all interno del dominio secondo una gri glia regolare equispaziata i valori di al variare del passo di griglia sono riportati in tabella 4 10 e in scala logaritmica in figura 4 25 parametro valore unit di misura L 8 m lato D 2 m E le3 GPa 1 1 97E 006 v 0 25 0 5 4 04E 007 P 2 N 0 25 6 59E 008 0 125 9 06E 009 Tabella 4 9 Parametri descrit Tabella 4 10 Valore dell errore tivi il problema della trave di Ti medio 6 per i punti interni al do moshenko usato per il confronto minio al variare del passo della numerico griglia di punti 129 4 L approccio meshless con il Metodo delle Celle
99. e esterno XOY come illustrato in figura 4 3 S3 S X S4 4 S5 2 S m Se gt x Figura 4 3 Distribuzione di nodi nel dominio in esame e definizione di una mesh locale attorno al generico nodo interno Pp con sei nodi satellite 5 Note le condizioni al contorno temperature imposte sorgenti di calore con centrate flussi di calore sul contorno ci si pone l obiettivo di determinare i valori della temperatura in ogni nodo L algoritmo prevede la definizione di una mesh di celle primali localmente nodo a nodo preso in esame un generico nodo polo P e si individuano i pi vicini nodi circostanti e si ordinano i o questi nodi in senso ad esempio antiorario rispetto al nodo polo in modo tale da avere una lista ordinata di N nodi satellite S 4l individuazione dei nodi circostanti pu avvenire con un criterio empirico di individuazione dei nodi contenuti all interno di una circonferenza di raggio R valutato a partire da informazioni circa la densit globale o locale dei nodi 83 4 L approccio meshless con il Metodo delle Celle e si definiscono N celle locali c le prime N 1 celle avranno vertici P S 5 41 l ultima avr vertici P Sw S1 e si provvede alla definizione della regione tributaria del nodo polo un modo immediato ad esempio quello di creare il complesso di celle duale locale a partire dal complesso primale locale appena definito 2 S1 S3 S S4 qd S5 Se
100. e generalizzata di una grandezza locale infatti g g z y z t Le grandezze globali e dunque nel mondo della realt discreta possono essere invece associate agli istanti di tempo ovvero essere riferite ad intervalli di tem po istanti ed intervalli sono le entit temporali di riferimento la cui definizione intuitiva perch corrispondente alla quotidiana esperienza di ciascuno l istan te quel riferimento temporale usato per collocare un evento ad esempio l arrivo dell autobus ad una fermata l intervallo tempo intercorrente tra due istanti invece il riferimento naturale per fenomeni che non sono puntuali ma che han no uno svolgimento nel tempo ad esempio il tempo necessario all autobus per raggiungere la fermata successiva Per fare un esempio di ambito fisico il vettore posizione di una particella r una variabile associata agli istanti mentre lo spostamento s una grandezza attribuita agli intervalli Al pari degli elementi geometrici anche istanti ed intervalli possono ricevere tanto un orientazione interna quanto un orientazione esterna 2 3 1 Orientazione interna La naturale orientazione interna di un intervallo di tempo T costituita dal verso assegnato all asse dei tempi un intervallo di tempo cio orientato dall istante precedente a quello successivo L orientazione degli istanti in qualche modo omologa all orientazione interna di punti e il suo significato convenzionale di scende
101. e locali pena l errato computo dei termini di sorgente 4 2 1 1 Nodi di bordo La procedura da seguire per i nodi di bordo rimane sostanzialmente inalterata come peraltro gia nell approccio tradizionale nel bilancio al nodo polo andranno tenuti presente anche eventuali flussi di calore applicati al contorno del dominio nella porzione p che interessa la cella duale del polo L equazione 4 3 assume dunque la forma b9 Qp Gp p 4 5 cET P 4 2 1 2 Alcuni test numerici Per verificarne l accuratezza l approccio descritto stato applicato alla soluzione della equazione di Laplace V w 0 su un domino piano quadrato di lato uni tario centrato nell origine di un sistema di riferimento cartesiano con assi paralleli ai lati nel dominio sono stati posizionati dei punti secondo una griglia regolare parallela ai lati Da questa griglia con spostamenti random e facendo in modo di avere un punto in ogni vertice e un certo numero di punti posizionati lungo cia scun lato si sono poi ottenute delle distribuzioni casuali di punti nel dominio Sui punti del contorno sono stati imposti i valori assunti dalla funzione armonica wr x y e cosy 4 6 soluzione del problema il cui grafico raffigurati in figura 4 7 Come parametro di valutazione della precisione della soluzione numerica viene Sper i punti posizionati sui lati lo spostamento avviene lungo il lato stesso mentre i punti posizionati sui vertici non subiscono al
102. e per la costruzione di modelli numerici basa ti sulle immagini diagnostiche in grado superare i maggiori limiti della procedura basata sul FEM l intento quello di individuare valutandone l appicabilit un processo di analisi maggiormente automatico di accuratezza compatibile con le esigenze cliniche e completabile in tempi accettabili Il lavoro qui esposto si colloca nel contesto dell attivit di ricerca del Laborato rio di Tecnologia Medica dell Istituto Ortopedico Rizzoli di Bologna struttura che ha finanziato il Dottorato e presso cui lo studio stato condotto Parte consi stente dell attivit del Laboratorio focalizzata sulla caratterizzazione dei tessuti ossei a varia scala e sull affinamento di procedure validate di simulazione nume rica per lo studio del comportamento meccanico dei segmenti ossei nella prospet tiva del loro trasferimento alla pratica clinica ortopedica Oggetto della ricerca stato dunque nello specifico l indagine attorno a metodi numerici innovativi al ternativi al FEM relativamente agli specifici problemi di meccanica dell osso Il focus applicativo uno degli elementi caratterizzanti il lavoro guidando sia la scelta delle alternative studiate che il processo di valutazione dell applicabilita e di accertamento dell accuratezza L elaborato esplora alcune ipotesi fondate su implementazioni numeriche inno vative derivate del Metodo delle Celle Il Metodo delle Celle CM un metodo
103. e spostamento secondo la direzione x e y rispettivamente del sistema di riferimento globale e con il pedice num un valore di soluzione numerica e con il pedice un valore di soluzione teorica Il calcolo stato condotto distribuendo dei punti all interno del dominio secondo una griglia regolare equispaziata i valori di 6 al variare del passo di griglia sono riportati in tabella 4 6 e in scala logaritmica in figura 4 18 passo di errore medio parametro valore unit di misura griglia 4satelliti 8 satelliti L 8 m 1 5 88E 004 D 2 m 0 5 5 3751E 004 4 78E 004 E le3 GPa 0 25 4 865E 004 4 34E 004 v 0 25 0 125 4 7666E 004 4 24E 004 P 2 N Tabella 4 6 Valore dell errore Tabella 4 5 Parametri descrittivi il medio per i punti interni al do problema della trave di Timoshen minio al variare del passo della ko usato per il confronto numerico griglia di punti 110 4 4 Interpolazione locale con funzioni di base radiale J 4 nodi 3 24 n 8 nodi 3 26 4 e 3 28 4 e BD 334 2 e 3 32 4 n 3 34 4 3 36 7 n E 3 38 T T T T T T T 1 1 0 9 0 8 0 7 0 6 0 5 0 4 0 3 0 2 log lato Figura 4 18 Test della trave di Timoshenko errore medio rispetto al passo di gri glia in scala logaritmica con approccio CM meshless a interpolazione polinomiale A scopo esemplificativo della generalit dell approccio si riportano in figu ra 4 19 e 4 2
104. egolare equisispaziata 1 5 Figura 4 20 Test della trave di Timoshenko campo di spostamenti otte nuto su una distribuzione casuale dall approccio meshless CM ad inter polazione polinomiale insieme di punti si definisce una interpolazione polinomiale dei valori nodali del la variabile di configurazione del campo studiato la regione tributaria del nodo in esame viene definita attorno a ciascun nodo e per semplicit di espressione formale dell equazione di bilancio viene assunta con una forma geometrica re golare e semplice le grandezze coinvolte nell equazione di bilancio sono ricavate tramite le leggi costitutive e le equazioni di struttura a partire da una funzione che interpola i valori nodali della variabili di configurazione Il vantaggio di questo approccio indubbiamente quello di non richiedere alcu na definizione nemmeno a livello locale di connettivit tra i punti non occorre ordinare i nodi per costruire un set di celle che condividano tutte il nodo in esa me n necessario alcun controllo circa la qualit delle celle locali cos create 112 4 4 Interpolazione locale con funzioni di base radiale per escludere eventuali celle distorte Un ultimo vantaggio rispetto alla formu lazione meshless basata sulla creazione di una mesh locale ha a che fare con la copertura del dominio ottenuta con l insieme delle regioni tributarie locali in presenza di sorgenti distribuite sorgenti di calor
105. ella prova sperimentale sono stati im pediti gli spostamenti di tutti i nodi superficiali corrispondenti alle regioni del segmento osseo immerse nel blocco di PMMA II carico invece stato applicato al nodo di vertice pi prossimo alla posizione registrata del punto di applicazione del carico in ciascuna configurazione di sollecitazione Figura 5 10 Sezione coronale a sinistra e trasversale a destra del modello ad elementi finiti 5 2 4 Il modello MCM Il modello MCM stato generato utilizzando come nodi i punti della griglia della scansione CT senza alcuna operazione di segmentazione La costruzione del modello stata diretta selezionata la regione di interesse con quota a partire dal piano di imposta della base di PMMA l unica fase di prepro cessing stata un sottocampionamento subsampling su una griglia regolare di passo 2 mm riducendo i numero di punti in ragione delle limitazioni di memo ria del calcolatore disponibile per condurre l analisi Ai punti a quota minima corrispondente al piano di imposta della base di PMMA del femore sono stati imposti spostamenti nulli replicando le condizioni di vincolo delle prove di cari co sperimentali Il carico corrispondente a ciascuna condizione di sollecitazione stato applicato al nodo del dataset CT pi prossimo alla posizione registrata del 143 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle punto di applicaz
106. elle Celle 2 2 2 4 Orientazione esterna di un punto Conservando l idea di una definizione convenzionale che permetta di conservare le propriet formali anche il punto pu essere dotato di un orientazione esterna Ipotizziamo allora di considerare un cubo dotato di orientazione interna e un punto posto al centro del cubo i sensi di rotazione definiti sulle facce del volume possono vedersi come sensi di rotazione delle semirette con origine nel punto ragionevole considerare orientazione esterna di un punto quella costituita dai dei sensi di rotazione delle semirette con origine nel punto gv i G orientazione interna Orientazione interna di un punto un punto orientato positivamente se un pozzo Orientazione interna di una linea la nozione base usata per dare un significato alla orientazione di tutti gli altri elementi geometrici Orientazione interna di una superficie un oriertazione compatibile dei suoi spigoli per es il verso per percorrere i suoi bordi Orientazione interna di un volume un orientazione compatibile delle sue facce E equivalente alla regola orientazione esterna Orientazione esterna di un volume scelta di un orientazione interna od esterna delle normali Orientazione esterna di una superf scelta di un orientazione interna o esterna delle normali Orientazione esterna di una linea l orientazione interna della superficie che attraversa la
107. elli di automazio ne partendo dai dati CT esse permettono la definizione di modelli che riprodu cono le peculiari caratteristiche morfologiche del segmento osseo studiato e che riportano la specifica distribuzione delle propriet meccaniche dei tessuti presen ti Non tutte queste metodologie alcune delle quali gi usate con successo in vari contesti sembrano per rispondere a tutti i requisiti di automazione di ge neralita di esattezza e robustezza richiesti particolarmente da un applicazione 9 1 La creazione di modelli numerici subject specific di segmenti ossei clinica Qualunque strada si adotti la costruzione di un modello subject secific deve comunque affrontare due fasi fondamentali 1 la generazione della mesh 2 la definizione delle propriet dei materiali Entrambe queste fasi si basano sull interpretazione e la manipolazione delle in formazioni rese disponibili dalle immagini CT 1 3 1 La Tomografia Computerizzata CT La Tomografia Computerizzata o CT Computed Tomography una metodolo gia generale che permette di ottenere immagini tomografiche ovvero di sezioni trasversali di un oggetto Questo possibile con una tecnica nota come retropro iezione con la quale tramite un computer si riesce a ricostruire l immagine di una sezione sulla base dei proiezioni fatte da angoli diversi Solitamente l espres sione Tomografia Computerizzata fa riferimento a immagini costruite su proiezioni p
108. ema consideriamo a titolo di esempio una costellazione di 4 satelliti 15i problema geometricamente visualizzabile con la disposizione di una superficie piana su tre appoggi perch l appoggio non sia instabile i tre punti devono essere distinti e non allineati 16 cfr 65 102 4 3 Interpolazione locale con funzioni polinomiali y y S2 S2 S1 O gt Si X Sa Pol x S3 P 00 aa I s Pad N e gt S3 S4 S4 Figura 4 15 Costellazione di 4 sa Figura 4 16 Costellazione di 4 sa telliti posti a uguale distanza dal telliti posti a uguale distanza dal nodo polo e allineati a coppie se nodo polo e allineati a coppie se condo direzioni parallele agli assi condo direzioni tra loro ortogonali del sistema di riferimento locale e inclinate di 7 4 rispetto agli assi del sistema di riferimento locale disposti come in figura 4 15 nel sistema cartesiano locale centrato nel polo la matrice Mc assume la forma Ea 400 1d 0 0 Me 1 0 d 0 4 39 1 d 0 d 0 1 0 d 0 ruotando questa costellazione di punti di 7 4 si ottiene la disposizione di figura 4 16 nel riferimento cartesiano locale la matrice dei coefficienti diviene ora 1 0 0 00 1 vd vd d Me 1 vd vd d 4 40 1 vd vd d 1 vd vd d In questo secondo caso il determinante della matrice Me risulta nullo questa seconda disposizione di punti risulta quindi inadatta Dal punto di vista geo metrico le due costellazioni sono equivalenti risultand
109. ente inverso sono identificate dalle stesse lettere cui sovrapposto il simbolo tilde Poich le caselle sono destinate ad ospitare una o pi grandezze coinvolte nella formulazione ne deriva che le variabili di configurazione si collocheranno nelle caselle della colonna di sinistra mentre le variabili di sorgente troveranno collo cazione nelle caselle della colonna di destra Le equazioni di struttura che legano 42 2 6 Impostazione del problema con il Metodo delle Celle variabili di configurazione variabili di sorgente equazioni costitutive Figura 2 6 Schema generale di un diagramma di Tonti per l illustrazione del le variabili delle entit geometriche n celle di riferimento e delle equazioni che definiscono le relazioni tra le grandezze coinvolte variabili dello stesso tipo sono rappresentate da frecce verticali le equazioni co stitutive che costituiscono invece una sorta di ponte tra le due categorie di grandezze rappresentate da frecce orizzontali od oblique che collegano caselle appartenenti a colonne opposte L equazione risolvente si ottiene componendo le equazioni di struttura e quelle costitutive ottenendo un legame fra variabile di configurazione e quella di sorgente del campo Lo schema proposto certamente uno schema generale e non detto che tutti i legami che possibile indicare debbano necessariamente essere sempre presenti Nel caso di problemi stazionari ad esempio il
110. ente la creazione della mesh parte dalla superficie di confine del tessuto osseo in forma tassellizzata Da questa base di partenza sono ancora possibili due metodi alternativi e rappresentare la superficie tassellizzata con una serie di superfici matema ticamente definite che nel loro complesso definiscano un solido topologica mente corretto Questa superficie pu allora essere discretizzata e la mesh ottenuta serve di appoggio per la creazione della mesh del volume Per descrivere la superficie tasselizzata si possono usare degli elementi di su perfic NURBS Non Uniform Rational B Splines generabili con sofisticati programmi di reverse engineering e sulle quali la creazione di una mesh di superficie risulta relativamente facile Importate le NURBS in un program ma FEM si provvede dapprima alla creazione della mesh sulle NURBS che funge da appoggio alla discretizzazione del volume con elementi tetraedri ci e la superficie tassellizzata ottenuta con la segmentazione pu essere utilizza ta senza ulteriori elaborazioni La metodologia chiamata grid based ef fettua sostanzialmente una proiezione della griglia di punti del dataset sulla superficie di contorno L assenza di codici commerciali che implementino questa procedura ne ha di fatto determinato la scarsit di utilizzo 17 1 La creazione di modelli numerici subject specific di segmenti ossei 1 3 3 Definizione delle propriet meccaniche Una volta completata la
111. eory source variables geometrical variables statical variables primal complex inner orientation dual complex outer orientation dimensions L dimensions MLT Ix P u displacement F volume force h relative displacement T surface force H displacement gradient A area vector strain tensor o symmetric stress tensor bulk strain p mechanical pressure e strain deviatoric o stress deviatoric b Burgers vector p stress vector potential H HT THT Poe Wo gt in Dt 93 3 Figura 2 8 Diagramma di Tonti per la formulazione finita dell elestostatica 45 2 Il Metodo delle Celle 46 Capitolo 3 Formulazione quadratica con il Metodo delle Celle Rispetto al Metodo degli Elementi Finiti il Metodo delle Celle utilizza un ap proccio metodologico completamente diverso l uso delle variabili di dominio in luogo delle funzioni di punto l attribuzione delle grandezze agli elementi geo metrici secondo i processi sperimentali utilizzati per la loro misura sono due elementi fondamentali che portano ad una formulazione discreta diretta senza dover fare ricorso alla formulazione differenziale stato dimostrato 30 che pur usando un approccio completamente diverso il CM usando i simplessi in grado di ottenere la stessa matrice di rigidezza costruita con il Metodo degli Ele menti Finiti Relativamente all uso dei simplessi la differenza tra i due metodi in termini di risultati ottenib
112. equazione di bilancio pu essere espressa come Que Oe Ore 2 1 dove Qrre la produzione netta della grandezza Q entro il volume V nell intervallo di tempo T Q l accumulo netto della grandezza Q entro il volume V nell intervallo di tempo T Q 5 il flusso uscente netto della grandezza Q attraverso la superficie di contor no del volume V nell intervallo di tempo T A sua volta la quantit accumulata Q costituita dalla differenza Or Qfinale Cade 2 2 fra la quantita contenuta entro V all istante finale e all istante iniziale dell inter vallo di tempo T su cui si considera il bilancio Quando il termine di produzione di Q sia nullo si parla di equazione di conser vazione sono equazioni di conservazione l equazione di bilancio della massa per un sistema chiuso e l equazione di bilancio dell energia di un sistema isolato Un altra specializzazione dell espressione generale dell equazione di bilancio si 40 2 6 Impostazione del problema con il Metodo delle Celle ha nel caso di sistemi stazionari ovvero sistemi che non hanno evoluzione nel tempo E il caso del bilancio statico delle forze che descrive l equilibrio di un corpo per il quale si ha un espressione generale del tipo F 0V F V 0 2 3 dove F V descrive le forze di volume agenti sul volume V mentre il termine F OV le forse di superficie che agiscono sul suo contorno Va notato che la forza pu essere vista co
113. ervallo duale pr SEB RRB LL da E istanti duali nici I I I istanti primali intervalli primali T T Figura 2 4 Orientazione esterna degli istanti e degli intervalli di tempo Si noti che sull asse dei tempi si sono individuate due classi entit gli istanti tem porali I tra cui sono individuati gli intervalli T costituiscono un primo ordine di elementi denominati primali Da questi elementi si sono ricavati gli istanti ed i corrispondenti intervalli T queste entit sono denominate duali poich corri spondenti alle entit primali e sono ad esse strettamente correlate Istanti e intervalli primali sono dotati di orientazione interna istanti e intervalli duali sono invece dotati di orientazione interna indotta dall orientazione delle corrispondenti entit primali 34 2 4 Attribuzione delle variabili agli elementi geometrici 2 4 Attribuzione delle variabili agli elementi geome trici Le grandezze fisiche globali sono associate ad elementi geometrici orientati In termini generali si pu affermare che la definizione stessa di una variabile ov vero le sue modalit di misura sperimentale forniscono l indicazione dell ente geometrico cui la grandezza riferita Circa invece l orientazione di questo ente geometrico interna o esterna si osserva che per molte variabili il cambiamento di orientazione dell ente geometrico comporta un cambiamento di segno della variabile stessa Queste affermazioni
114. ess con il Metodo delle Celle FEM risulta inadeguato in caso di grandi deformazioni con la perdita di continui t del dominio Lo sforzo tangibile di ricorrere ad approcci meshless per raggiungere simula zioni complesse difficoltosamente approntabili altrimenti ravvisabile in molti lavori riportati in letteratura Wong et al 53 presentano ad esempio un am biente per l applicazione di approssimazioni definite element free nell ambito della cardiologia computazionale in cui viene presentata l implementazione di un modello meshless di un modello fisiologico del cuore costruito da immagini medicali Il lavoro riporta delle comparazioni con il FEM su geometrie semplifi cate e dimostra la plausibilit fisiologica di un modello di cuore umano costruito da immagini di risonanza magnetica Chun et al 54 realizzano un modello MGED per lo studio della crescita della placca interna dei vasi sanguigni il modello basato su tre immagini di risonanza magnetica distanziate nel tempo bidimensionale e adotta un modello semplifi cato sia per le propriet meccaniche del vaso che per il modello di accrescimento del deposito Horton et al 55 affrontano la simulazione della deformazione di tessuti mol li in tempo reale con una formulazione lagrangiana meshless basata sulle diffe renze centrali con integrazione esplicita per nonlinearit di carattere geometrico il focus dell applicazione la brain indentation problema rispett
115. essivo prevedeva l applicazione sulla testa femo rale di una forza pari a circa il peso corporeo 800 N per evitare condizioni di 2polimetilmetacrilato 139 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle Figura 5 8 Rappresentazione schematica del campione montato sulla macchina di prova in posizione di flessione a 18 e rappresentazione schematica delle confi gurazioni di carico il carico applicato verticalmente neutral LC1 e con diverse inclinazioni rispetto all asse biomeccanico nel piano frontale adduzione 3 LC2 e 24 LC4 e nel piano sagittale 3 estensione LC3 e 18 flessione LC5 danno meccanico nel campione secondo 5 diverse configurazioni di carico cor rispondenti alle direzioni estreme raggiunte dalla forza di reazione all anca nello svolgimento di normali azioni motorie 79 La prova stata condotta in controllo di carico consentendo cos lo sviluppo di eventuali effetti di creep il carico mas simo raggiunto quasi staticamente con incrementi di 100 N stato mantenuto per circa 30 secondi Sono state eseguite cinque ripetizioni di misura ad interval li di almeno 10 minuti per consentire il recupero di eventuali effetti viscoelastici indotti La procedura descritta ha consentito di misurare con ottima ripetibilit le deformazioni controllare la linearit della risposta all applicazione del carico quantificare l entit di fenomeni
116. ette vs sperimentali per le varie condizioni di carico 6 4 Comparazione sistematica con modello FEM di ri ferimento Al fine di approfondire il comportamento del modello MCM stato disposto un confronto puntuale tra un modello MCM riproducente uno dei femori in una del le configurazioni di carico e il corrispondente modello FEM Il confronto indaga il comportamento relativo tra i due modelli confronto nel quale il FEM assunto come riferimento non gi quale soluzione esatta quan to piuttosto come migliore strumento validato attualmente disponibile nel ripro durre numericamente le misure sperimentali Va anticipato che la valutazione del comportamento relativo del modello MCM rispetto al FEM deve tenere presente 165 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle 600 9 pe MCM 400 M 200 5 He exp 600 400 200 F 200 400 600 5 d P adduzione 3 200 y 1 02 x 4 55 DA R 0 95 E RMSE 61 88 pe 10 81 aoo Max Err 179 56 pe 31 38 aor 600 800 600 400 200 pe MCM Sa ue exp 400 200 neutral y 1 03 x 7 90 R 0 94 RMSE 79 36 pe 11 09 Max Err 237 90 pe 33 24 pe MCM 200 400 600 800 200 400 adduzione 8 y 1 09 x 1 67 R 0 95 RMSE 62 54 pe 12 31 Max Err 191 62 pe 37 73 200 400 600 800 estensione 3 y 1 07 x 7 05 R 0 93 RMSE
117. finizione di un modello densita al modulo di Young le due operazioni non sono commutative pertan to una strada alternativa e non equivalente quella di passare in ogni punto del dataset dai valori CT ai valori del modulo di Young associando poi ad ogni elemento la media dei valori del modulo di Young ottenuti nei voxel componenti Avendo generato una mesh strutturata possibile realizzare un accurata map patura delle proprieta dei diversi tessuti dato che le facce degli elementi possono seguire debitamente i confini fra i vari tessuti ossei In questo caso la densita media in un singolo elemento finito definita sulla base dei voxel che cadono al proprio interno La procedura pu avere diversi livelli di automazione ma in ogni si basa implicitamente sulla conoscenza a priori della topologia della mesh 15 13 16 Quando si sia invece costruita una mesh non strutturata le facce degli ele menti non sono allineate n con gli assi del dataset CT n con le superfici del tessuto osseo poich improponibile procedere ad una mappatura manuale del le propriet meccaniche sugli elementi della mesh necessario escogitare dei metodi pi sofisticati Se la mesh generata partendo dagli stessi dati la mesh e la distribuzione delle densit sono perfettamente registrate nello spazio L unico problema quello di stabilire come trasferire agli elementi del modello i valori di questa non omogenea distribuzione delle propriet
118. formazione Posizionati all interno della geometria del corpo un certo numero di punti nodi riferiti ad un riferimento cartesiano globale lo scopo determinare le componenti del vettore spostamento in ciascuno punto 106 4 3 Interpolazione locale con funzioni polinomiali La chiave del metodo la scrittura del bilancio delle forze ad ogni nodo definita una regione tributaria 2 associata ad ogni nodo i si impone l equivalenza tra le forze di volume Bo agenti su tale regione e le forze di superficie Typ agenti sul suo contorno 02 Bot YO Tani 0 4 45 L equazione vettoriale contempla le componenti delle varie forze riferite al si stema di riferimento cartesiano globale Considerato un generico nodo P interno al dominio si individui tra in nodi prossimi una costellazione di nodi satellite S per semplicit di illustrazione si assumer nel seguito che la costellazione sia composta da 5 satelliti si adotti inol tre un sistema di riferimento cartesiano locale con origine nel nodo polo e assi paralleli al sistema di riferimento globale 4 3 2 1 Interpolazione locale Nell intorno del nodo polo il campo degli spostamenti viene approssimato con 4 46 2 2 fa A 03 Asy 07X a9XY a11y 2 2 U x y az asx age asy aioty 012Y 4 3 2 2 Deformazioni Le componenti del tensore delle deformazioni vengono ricavate tramite le oppor tune derivate delle componen
119. funzione analiticamente definita ad esempio un polinomio di grado opportuno e derivare analiticamente le grandezze coinvolte nell equazione di bilancio attraverso le equazioni di struttura e costitutive applicate analiticamente alle funzioni interpolanti adottate e scrivere l equazione bilancio del nodo polo relativamente ad una regione tributaria di forma elementare circonferenza sfera quadrato cubo cen trata sul polo polo stesso Da quanto esposto resta lecito chiedersi quanto un simile approccio sia dal punto di vista della categorizzazione affiliabile al Metodo delle Celle del Meto do delle Celle viene infatti conservata solamente l idea di bilancio su una regione tributaria di nodo ad estensione discreta per il resto la formulazione algebrica delle equazioni coinvolte risulta appannata da un significativo uso di equazioni integro differenziali pure la filosofia del riferimento delle variabili agli elementi geometrici risulta sfumata A lungo si potrebbe quindi discutere sulla la propriet con cui la letteratura abbia accettato almeno relativamente alla prima delle due formulazioni che verran no di seguito proposte di annoverare questi approcci nella famiglia del Metodo delle Celle la questione che pure tutt altro che marginale appare in questo contesto di secondario rilievo oggetto della ricerca non infatti specificamente l indagine circa l applicabilit del Metodo delle Celle ai problemi di biomeccani ca il l
120. garitmica di errore griglia di punti e passo di griglia deformazioni superficiali e i valori predetti dal modello numerico nei punti cor rispondenti ha costituito la base di valutazione dell accuratezza dell intera pro cedura di modellizzazione numerica Questo esperimento ha dunque costituito una prima importante doppia verifica dell approccio MCM la valutazione di ac curatezza avvenuta per confronto con le misure sperimentali mentre i risultati ottenuti dal FEM a livello aggregato hanno costituito una base di valutazione comparativa rispetto alla procedura di simulazione di riferimento 5 2 1 L esperimento di riferimento Su un femore destro proveniente da un donatore maschio di 51 anni peso 75 kg altezza 175 cm sono state applicate 13 rosette estensimetriche triassiali posi zionate su ciascuno dei quattro aspetti anatomici anteriore posteriore mediale laterale su quattro diversi livelli in direzione prossimale distale testa collo metafisi diafisi Nella regione diafisaria sono stati applicati soltanto due esten simetri in posizione laterale e mediale ove si suppone si misurino deformazioni pi significative per limitare il numero di canali di acquisizione L applicazione di un estensimetro nell aspetto laterale del collo stata parimenti omessa a cau sa dell estrema irregolarit della superficie Si faccia riferimento alla figura 5 7 per un identificazione pi chiara del posizionamento dei sensori La superfic
121. generato a partire dal dataset CT seguendo una procedura che pu essere cos schematizzata e estrazione del modello virtuale della superficie tridimensionale dell osso e generazione della mesh di elementi finiti e mappatura delle propriet meccaniche modello costitutivo elastico isotro po disomogeneo del tessuto osseo sulla mesh Per quanto riguarda il primo passo i contorni esterni della superficie del campio ne sono stati estratti dalle sezioni tomografiche a mezzo di un software semiau tomatico basato su un algoritmo di identificazione dei bordi 80 I contorni cos ottenuti sono stati interpolati per ottenere la superficie esterna del femore usan do una triangolazione di Delaunay 5 La superficie tridimensionale tassellizzata cos ottenuta stata convertita in un modello NURBS Non Uniform Rational B Spline con il software di reverse engineering Geomagic Studio v 6 Raindrop Geomagic Inc USA Il secondo passo di modellazione stato realizzato generando in maniera au tomatica con il pacchetto software HyperMesh Altair Engineering Inc USA che implementa un algoritmo advancing front una mesh di tetraedri a dieci nodi del femore stata posta attenzione a che sulla superficie del modello fosse generato uno strato regolare di tetraedri in modo da minimizzare l inevitabile effetto di volume parziale che tende a sottostimare la densit e quindi il modulo elastico assegnata al modello nelle regioni di confine e pa
122. gions 3 In entrambi i casi i valori ottenuti non possono essere considerati privi di errore Successivamente i valori di densit del tessuto osseo sono correlati con 22 1 3 La definizione di un modello i corrispondenti valori di modulo di Young utilizzando equazioni di derivazio ne empirica ottenuti da indagini sperimentali I coefficienti di queste equazioni sono affetti da un incertezza dovuta a una significativa dispersione di queste mi sure sperimentali 26 27 28 29 Le inevitabili incertezze legate ai vari passi di creazione del modello affliggo no sicuramente l accuratezza dei valori ottenuti con il calcolo Sfortunatamente questi modelli subject specific non possono essere convalidati da misure dirette poich le tecnologie attualmente disponibili non permettono misure di tensione o di deformazione in vivo che non siano invasive La valutazione dell accuratezza del modello pu essere ricavata con due ope razioni separate dapprima si pu accertare l accuratezza intrinseca dei metodi messi a punto rispetto ai risultati sperimentali di test in vitro poi pu si eseguire un analisi di sensibilit per capire quale sia l effetto delle incertezze delle misure rilevate in vivo Mentre la validazione in vitro dei modelli ad elementi finiti di segmenti ossei stata realizzata da molti autori ancora sconosciuta l influenza che gli errori associati alle misure prese in vivo possano avere sui risultati prodotti da u
123. goli e punti a cui poter attribuire le grandezze da collocare nello spazio Questa strut tura costituita dai complessi di celle In termini pratici l operazione non nuova si tratta di dividere il dominio in tan te porzioni celle di forma arbitraria e di dimensioni opportune Questo primo complesso di celle detto complesso primale esso sin qui analogo all insieme di elementi delle mesh che si usano comunemente quando si adopera un meto do numerico tradizionale che discretizza le equazioni differenziali Nonostante questa iniziale analogia si preferisce utilizzare il termine complesso di celle in luogo di mesh in quanto tutti gli elementi geometrici che formano il complesso sono coinvolti nella formulazione discreta Vi poi un altra differenza fra mesh tradizionalmente intesa e i complessi di celle dato un complesso di celle di forma qualunque si pu individuare in ma niera arbitraria un punto all interno di ogni cella Congiungendo tali punti interni per ogni coppia di celle adiacenti si ottiene un secondo complesso di celle A que sto secondo complesso di celle si d il nome di complesso duale Vi una stretta corrispondenza reciproca dualit tra i due complessi di celle tale da giustificare l attributo di duale nella denominazione del secondo siste ma di celle innanzitutto ad ogni elemento del complesso primale corrisponde un opportuno elemento del complesso duale secondo il legame illustrato ne
124. i Enrico Schileo Benedikt Helgason Luca Cristofolini and Marco Viceconti The material mapping strategy influences the accuracy of ct based finite element models of bones An evaluation against experimen tal measurements Medical Engineering amp Physics 29 9 973 979 November 2007 Willi A Kalender Dieter Felsenberg Harry K Genant Manfred Fischer Jan Dequeker and Jonathan Reeve The european spine phantom a tool for standardization and quality control in spinal bone mineral measurements by dxa and qct European Journal of Radiology 20 2 83 92 1995 Elise F Morgan Harun H Bayraktar and Tony M Keaveny Trabecular bone modulus density relationships depend on anatomic site Journal of Biomechanics 36 7 897 904 2003 Dieter Christian Wirtz Norbert Schiffers Thomas Pandorf Klaus Raderma cher Dieter Weichert and Raimund Forst Critical evaluation of known bo ne material properties to realize anisotropic fe simulation of the proximal femur Journal of Biomechanics 33 10 1325 1330 2000 Marco Viceconti Sigbjorn Olsen Lutz P Nolte and Kim Burton Ex tracting clinically relevant data from finite element simulations Clinical Biomechanics 20 5 451 454 June 2005 Enrico Schileo Fulvia Taddei Andrea Malandrino Luca Cristofolini and Marco Viceconti Subject specific finite element models can accurately pre dict strain levels in long bones Journal of Biomechanics 40 13 2982 2989 2007 Enri
125. i z iui a ie Zi B o iM Tj Yi AF Yi Bi Zj t UM vi Yj di B Ds i Yj Tj Yi Zk Le Yi Ti Yu Zj Lj Ye Tk Yj zi B essendo A E xi eae 21 zi p i Um yi z 2 uz yi 2 z ys Ys zr 2i Yk Ys z zi x yi yi zk zi B xj i ye i xx Y Er Tj ul z x j Y T t ee xi ut xj ca w Spe Yj j T1 Yi Zk i T1 Yk T1 Tk Yil zj i Tj Yk Tk T1 Ys zi 50 3 5 3 6 3 7 3 8 3 9 3 10 3 11 3 12 3 13 3 14 3 15 3 16 3 17 3 18 3 1 Problemi scalari tridimensionali 3 1 1 2 La funzione interpolante Nel sistema di riferimento delle coordinate affini locali il polinomio che interpola i valori delle temperature nodali all interno della cella assume la forma Teac Stl Seg nie a ah Qi Qj ak 3 19 Am Un Qo ap I dieci coefficienti incogniti si determinano imponendo che il polinomio assuma ai nodi i valori noti _ or 0 1 2 1 2 0 0 1 2 0 H__ ji jd eC ee ee Li 0 0 0 0 1 0 0 1 0 0 1 2 0 1 0 0 0 1 2 0 1 4 1 2 1 2 0 0 0 0 0 1 4 0 0 0 0 OO SS Sa OS Si 1 4 a OD 5 O O O SL KH 4 0 0 0 0 1 0 so ah Th Oe 1 0 i Ti 0 0 1 j 1 Qk Th 1 4 0 0 o ia o 0 0 1 4 ki 1 4 0 1 4 7 0 1 4
126. i deformazione ai centroidi del modello FEM a sinistra rappresentato il diagramma cumulativo della frequen za dell errore percentuale delle predizioni MCM rispetto ai valori FEM a destra tiportata la distribuzione delle differenza con segno Il grafico di figura 6 13 rappresenta il diagramma cumulativo della frequen za di si evince una sostanziale concordanza tra i due modelli espressa dal fatto che circa il 60 dei valori ha uno scarto percentuale dai valori FEM entro il 10 entro uno scarto del 50 compreso il 92 30 dei valori confrontati Il primo di questi due dati assume particolare rilievo se si considera che lo scarto rispetto al FEM dell ordine dell errore che caratterizza le predizioni FEM ri spetto alle misure sperimentali In altre parole entro questo scarto non dato di sapere quale tra i due modelli FEM e MCM sia il pi coerente con la realt sperimentale Le differenze nei valori che caratterizzano 1 8 dei confronti possono essere ra gionevolmente imputate a tre fattori principali x e un primo elemento di disturbo rappresentato dall errore di posiziona mento che caratterizza il modello FEM rispetto al riferimento del dataset CT questo fatto origina in fase di segmentazione quando viene definita la superficie tassellizzata del segmento osseo poi oggetto di descrizione ma tematica tramite patch di superfici NURBS Al termine di questa procedura il posizionamento del modello a meno di
127. i esterne di un segmento consi derato nel piano il segmento orientato esternamente da un verso di attraversa mento nello spazio l orientazione esterna viene definita da un senso di rotazione attorno ad esso in uno spazio unidimensionale l orientazione invece costituita da due frecce che ne stabiliscono idealmente come positiva la trazione ovvero la compressione Figura 2 2 L orientazione esterna di una linea dipende dalle dimensioni dello spazio in cui la si considera 2 2 2 2 Orientazione esterna di una superficie La pi nota orientazione esterna di una superficie quella che viene fatta nello spazio attribuendole un verso di attraversamento in relazione a questo tipo di orientazione che assume significato il flusso di una grandezza attraverso una superficie In un ambiente bidimensionale un elemento di superficie riceve un orientazione esterna quando sia attribuita una direzione positiva di attraversamento del suo contorno ovvero quando il suo contorno abbia ricevuto un orientazione esterna 2 2 2 3 Orientazione esterna di un volume L orientazione esterna di un volume ristretta evidentemente al solo caso tridi mensionale non sussistendo volumi immersi nel piano costituita da un ver so di attraversamento della sua superficie di contorno Da notare che anche in questo caso l orientazione esterna di un volume viene indotta dall orientazione esterna della superficie di contorno 31 2 Il Metodo d
128. i nodi viene ritenuta un elemento imprescindibile tanto che Idelsohn et alii 47 e 48 affermano Un metodo meshless un algoritmo in cui 1 la definizione delle funzioni di forma dipende solamente dalla posizione dei nodi 2 la valutazione della connettivit locale ha complessit computazionale lineare ri spetto al numero dei nodi Questi requisiti giustificano quanti affermano che in realt sotto il profilo della efficienza nella definizione della connettivit non tutti i metodi siano real mente meshless Lo sviluppo degli algoritmi comunque ancora in corso e come alle origini i limiti del FEM oggi sono le difficolt e le limitazioni degli algoritmi meshless storici ad essere l ostacolo che le nuove formulazioni si impegnano a supera re come ben evidenziato in 49 Allo stato attuale gli approcci piiu promettenti paiono essere rappresentati dai metodi MFEM XFEM e GFEM 4 1 1 Utilizzo di approcci meshless in ambito biomeccanico Nello specifico dell ambito biomeccanico computazionale l applicazione di ap procci meshless resta tutt ora limitata a pochi studi malgrado sia gli aspetti mor fologici delle geometrie che la complessit dei problemi trattati rendano a priori interessante il loro utilizzo I principali ambiti di applicazione sono la simulazione chirurgica 50 51 e la simulazione del cedimento dei tessuti molli 52 ambiti entrambi in cui il lofr 48 77 4 L approccio meshl
129. i punti corrispondenti ai siti di misura sperimentale Le predizioni di un modello FEM replicante la stes sa condizione sperimentale fornisce anche elementi di valutazione comparativa rispetto al gold standard numerico di riferimento Il capitolo 6 infine dedicato alla validazione della formulazione MCM la cui applicabilit stata verificata nel capitolo precedente Come riferimento viene adottata una campagna di misure sperimentali condotte su 8 differenti campio ni di femore in 6 distinte configurazioni di carico le misure sperimentali di de formazione in 15 differenti siti anatomici costituiscono l elemento di verifica di precisione nell identificazione della realt sperimentale Definita una metrica di accuratezza globale del modello le prestazioni registrate dalle analisi MCM so no valutate anche per comparazione con i corrispondenti parametri ottenuti con i modelli FEM replicanti i medesimi esperimenti Il confronto puntuale e sistemati co con i valori di deformazione principale negli elementi di uno dei modelli FEM offrono infine uno strumento di comprensione del comportamento del modello MCM rispetto a quello del metodo numerico validato di riferimento Capitolo 1 La creazione di modelli numerici subject specific di segmenti ossei La conoscenza delle tensioni che si producono all interno di un segmento osseo pu risultare di grandissima importanza per una molteplicit di applicazioni sia di carattere clinico che nell
130. ia di tipo lineare che di tipo esponenziale Si assume che una singola equa zione resti valida per tutto l intervallo di valori di densit L equazione adottata nella forma E atop 1 4 nella quale E il valore del modulo di Young assegnato uniformemente all inte ro elemento n della mesh p il corrispondente valore di densit e a b e c sono i coefficienti adottati dall utente Teoricamente il procedimento pu assegnare pro priet differenti ad ogni elemento della mesh siccome alcuni codici FEM com merciali hanno limitazioni sul numero massimo di materiali presenti il codice offre all utente la possibilit di ridurre il numero di valori per il modulo E impo stando un valore soglia per AZ con una procedura in grado di assicurare che ad 21 1 La creazione di modelli numerici subject specific di segmenti ossei ogni elemento non venga attribuito un valore inferiore a quello individuato con l espressione 1 4 1 3 4 Iproblemi dell approccio FEM mesh I procedimenti utilizzati per estrarre le informazioni geometriche e materiali dai dati CT sono spesso affette da errori non trascurabili questi errori si propagano in modo ignoto nei vari passi della creazione del modello e affliggono in modo difficilmente accertabile l accuratezza dei calcoli Errori di risoluzione dello scanner CT La prima fonte di errore e distorsione costituita dalla risoluzione a dalla qua lit delle immagini del dataset utilizzato carat
131. iabili energetiche sono ottenute dal prodotto di una variabile di confi gurazione per una variabile di sorgente Sono variabili energetiche il lavoro la potenza e l energia nelle sue diverse forme l energia potenziale l energia cinetica l energia interna Il Metodo delle Celle tratta variabili globali e trae grandi vantaggi teorici dal l ultima delle tre classificazioni proposte Lo scopo di una classificazione quel lo di ordinare gli elementi classificati chiarendo le relazioni che li accomunano e quelle che li differenziano Le classificazioni proposte sono di tipo funziona le e sono di grande aiuto nell impostazione della formulazione discreta delle equazioni di campo Questo fatto legato alla stretta relazione che si instaura tra le classi di variabili e orientazione degli elementi geometrici cui le variabili sono riferite Come verr illustrato nel seguito esiste in sostanza una partico lare corrispondenza tra le variabili fisiche ed attributi degli enti geometrici di appoggio Le variabili fisiche globali siano esse variabili di configurazione di sorgente o energetiche vengono riferite alle strutture geometriche elementari dello spazio punti linee superfici e volumi e del tempo istanti ed intervalli di tempo Il segno di queste variabili dipende dall orientazione attribuita a questi elementi spaziali e temporali di riferimento Considerando ad esempio un flusso attraverso una superficie del tutto eviden
132. icando la precisione della soluzione ottenibile in un problema di riferimento standard e comparando i risultati sia con la soluzione teorica del problema che con i valori ottenuti da un analisi FEM di confronto 3 1 Problemi scalari tridimensionali In questo paragrafo verr presentata l implementazione del Metodo delle Celle con interpolazione quadratica nella soluzione dell equazione di Poisson in domi ni tridimensionali Come problema di riferimento viene assunta la conduzione termica in regime sta zionario Dopo l illustrazione dell implementazione verr descritto il test di valu tazione dell accuratezza nella soluzione di un problema di controllo di cui nota la soluzione analitica comparando i risultati con quanto ottenuto da un analisi gemella con identiche mesh interpolazione e condizioni al contorno condotta con un codice commerciale FEM 3 1 1 Interpolazione quadratica Il campo delle temperature viene approssimato all interno di ogni cella con un polinomio di secondo grado T a y z a bx cy dz exy fyz gzz 4 ha iy 12 3 1 La funzione che ha 10 coefficienti richiede l uso di un tetraedro a 10 nodi la cella viene completamente definita dai suoi 4 vertici e dai punti medi dei suoi 6 lati si vuole interpolare la temperatura in termini dei valori ai nodi Al fine di poter disporre di una formulazione assolutamente generale ed indipendente dal la geometria della singola cella conveniente adotta
133. icina di cui fu insignito nell anno 1979 14 1 3 La definizione di un modello struttura ossea produce una pila di immagini che raffigurano sezioni trasversali del segmento studiato Noto il riferimento geometrico che permette di collocare le immagini nello spazio la pila di immagini riporta i valori di un campo scalare nei punti di una griglia spaziale regolare il campo scalare quello dei coefficienti di attenuazione definiti indipendentemente dai confini di separazione fra tessuti di organi differenti La soluzione agli elementi finiti delle equazioni differenzia li che descrivono un problema richiede per che il dominio di integrazione sia chiuso e ben definito In questa struttura di dati denominata dataset la distinzione fra i tessuti non determinabile a priori e dunque la definizione della geometria del segmento os seo essenziale ai fini di una generazione automatica di una mesh non pu essere immediata Allo stato attuale sono state definite due strategie fondamentali per l individua zione della geometria e dunque per poter provvedere alla generazione automa tica della mesh di un segmento osseo 1 una procedura di sogliatura generalizzata con cui si individuano i punti che potrebbero appartenere al tessuto osseo 2 la definizione di una superficie di confine tra tessuto osseo e tessuti circo stanti attraverso un processo denominato segmentazione 1 3 2 1 Il processo di sogliatura automatica
134. ie di applicazione di ciascuno dei sensori stata preparata secondo una procedura va lidata in letteratura 78 per le prove su campioni idratati di osso cadaverico Le 138 5 2 Il test esplorativo 13 rosette estensimetriche applicate erano del modello KFG 3 120 D1711L3M2S Kyowa Tokyo Japan aventi una lunghezza di griglia di 3 mm La tensione applicata stata limitata a 0 5 V per griglia per evitare effetti di surriscaldamen to Le deformazioni sono state acquisite ad una frequenza di 10 Hz con un filtro passa basso in frequenza fissato a 1H z AN 9 MNO AI gt Lig M1 M5 Le M Figura 5 7 Schema della disposizione degli estensimetri da sinistra vista ante riore laterale posteriore mediale La parte distale del femore stata inclusa in una base di afferraggio di PM MA per realizzare un vincolo ad incastro nel posizionamento sulla macchina di carico secondo gli schemi di seguito descritti La macchina di prova era equipag giata con una cella di carico di portata 5kN Mod 8502 Instron Canton MA USA Il campione stato montato sulla cella di carico per mezzo di una piatta forma che consente l inclinazione del campione secondo diversi angoli di prova Il carico veniva trasmesso dall attuatore attraverso slitte a basso attrito per evita re la trasmissione di indesiderate componenti orizzontali di carico Uno schema dell applicazione del carico riportato in figura 5 8 Lo schema di carico compl
135. ili Esistono oggi numerosi codici commerciali che integrano formulazioni FEM per la soluzione di problemi complessi in cui intervengano problemi afferenti a distinti problemi fisici Pur evoluto talmente da consentire analisi molto raffinate e complesse il FEM non privo di limiti che in taluni ambiti fanno s che si siano cercate strade alternative o parallele a maggiore flessibilit o di minor onere computazionale Una delle fasi pi onerose di una analisi FEM la discretizzazione delle geome trie in elementi con la creazione della cosiddetta mesh dividere una geometria in 73 4 L approccio meshless con il Metodo delle Celle elementi contigui riproducendone le caratteristiche geometriche con sufficiente approssimazione e disponendo di elementi regolari infatti un problema tutt al tro che banale specie relativamente ai casi tridimensionali l operazione ha un esito non sempre garantito e spesso risulta essere un punto critico ai fini del la precisione dei calcoli effettuati La generazione della mesh una operazione spesso computazionalmente onerosa e dispendiosa in termini di tempo fatto che assume rilievo in tutti quei problemi dalla fluidodinamica a superficie libera ai fenomeni di frattura dalle grandi deformazioni alle interazioni fluido struttura in cui il fenomeno studiato comporta una variazione della geometria obbligando ad ogni passo dell analisi alla rigenerazione della mesh sulla geometria aggior nata
136. ili risiede solo nella differente costruzione del termine noto del sistema lineare che descrive la soluzione approssimata del problema In particolare i due metodi sono caratterizzati dallo stesso valore di ordine di con vergenza Quando per si adotti la formulazione quadratica il Metodo delle Celle ha di mostrato di poter ottenere il quarto ordine di convergenza superiore al valore ordinariamente ottenibile dal FEM con la medesima interpolazione Questo possibile grazie ad una particolare costruzione del complesso di celle duale otte nuta appoggiando la divisione delle celle primali ai punti di Gauss dei suoi lati Questa propriet stata verificata relativamente alla soluzione numerica dell e quazione di Poisson per domini bidimensionali Resta l interrogativo circa il fatto che questa possibilit di ottenere il quarto ordine di convergenza permanga nella soluzione dei problemi tridimensionali In questo capitolo verr presentata l implementazione dell interpolazione qua dratica per la soluzione dell equazione di Poisson in domini tridimensionali con il Metodo delle Celle verr dimostrato come ancora un opportuna costruzione 47 3 Formulazione quadratica con il Metodo delle Celle del complesso duale permetta di ottenere l ordine di convergenza 4 In conclusio ne l uso dell interpolazione quadratica verr presentato nella sua estensione allo studio dei corpi elastici in regime di piccoli spostamenti verif
137. in Fluids 20 1081 1106 1995 38 J M Melenk and I Babuska The partition of unity finite element method Basic theory and applications Computer Methods in Applied Mechanics and Engineering 139 289 314 1996 39 E Onate S Idelsohn O C Zienkiewicz and R L Taylor A finite point method in computational mechanics applications to convective transport and fluid flow International Journal for Numerical Methods in Engineering 39 22 3839 3866 1996 40 T J Liszka C A M Duarte and W W Tworzydlo hp meshless cloud me thod Computer Methods in Applied Mechanics and Engineering 139 1 4 263 288 1996 181 41 N Sukumar B Moran and T Belytschko The natural element method International Journal for Numerical Methods in Engineering 43 839 887 1998 42 S Shen S N Atluri The Meshless Local Petrov Galerkin MLPG Method Tech Science Press 2002 43 T Strouboulis K Copps and I Babuscaronka The generalized finite ele ment method an example of its implementation and illustration of its performance International Journal for Numerical Methods in Engineering 47 8 1401 1417 2000 44 Sergio R Idelsohn Eugenio Onate Nestor Calvo and Facundo Del Pin The meshless finite element method International Journal for Numerical Methods in Engineering 58 6 893 912 2003 45 N Sukumar N Moes B Moran and T Belytschko Extended finite ele ment method for three dimensional crack modelling Intern
138. infine carattere applicativo indagando applicabilita e accuratezza di una delle ipotesi formulate Il capitolo 1 presenta il tema della definizione dei modelli numerici subject specific di segmenti ossei costruiti su dati di immagini diagnostiche vengono descritti gli strumenti disponibili e le metodologie ad oggi pi affidabili analiz zando le limitazioni e le incertezze connaturate sia ai metodi di manipolazione che alla fonte di informazione Nel capitolo 2 viene presentato il Metodo delle Celle conferendo particolare rilievo all architettura logica e formale su cui si fonda la filosofia della formu lazione finita cuore del Metodo delle Celle sinteticamente esposta con alcune digressioni di carattere tassonomico sulle possibili classificazioni delle variabili fisiche sull orientazione degli elementi geometrici e temporali e sui complessi di celle Nel capitolo 3 viene presentata l estensione della formulazione quadratica del Metodo delle Celle ai problemi tridimensionali la descrizione dell impostazione formale e implementativa relativamente ai problemi di riferimento trasmissione del calore in regime stazionario ed elasticit lineare completata da alcuni test su casi di riferimento in cui la soluzione numerica viene confrontata con quella teorica disponibile in forma chiusa Nel capitolo 4 vengono esplorate alcune ipotesi di approccio meshless del CM una breve disamina di metodi meshless in ambito differenziale viene c
139. ion with radial basis functions Advances in Engineering Software 38 5 320 327 May 2007 75 J Park and I W Sandberg Universal approximation using radial basis function networks Neural Computation 3 2 246 257 1991 76 Fulvia Taddei Martino Pani Luigino Zovatto Enzo Tonti and Marco Vi ceconti A new meshless approach for subject specific strain prediction in long bones Evaluation of accuracy Clinical Biomechanics 23 9 1192 1199 2008 77 Fulvia Taddei Luca Cristofolini Saulo Martelli H S Gill and Marco Vice conti Subject specific finite element models of long bones An in vitro eva luation of the overall accuracy Journal of Biomechanics 39 13 2457 2467 2006 78 M Viceconti A Toni and A Giunti Experimental Mechanics Technology Transfer Between High Tech Engineering and Biomechanics chapter Strain guage analysis of hard tissues factors influencing measurements Elsevier 1992 79 G Bergmann G Deuretzbacher M Heller F Graichen A Rohlmann J Strauss and G N Duda Hip contact forces and gait patterns from routine activities Journal of Biomechanics 34 7 859 871 2001 184 80 81 82 83 84 85 86 87 88 Debora Testi Cinzia Zannoni Angelo Cappello and Marco Viceconti Border tracing algorithm implementation for the femoral geometry recon struction Computer Methods and Programs in Biomedicine 65 3 175 182 2001 Fulvia Tadde
140. ione del carico sperimentale Attorno ad ogni nodo del dataset stata costruita la mesh locale di tetraedri costruiti con i nodi prossimi nelle direzioni di allineamento della griglia di punti A tutte le celle della mesh locale stato associato un valore unico di modulo di elasticit caratterizzante il legame costitutivo elastico lineare isotropo ad esse as sociato Il valore del modulo di elasticit stato determinato come media pesata dei valori determinati ai vertici tramite le relazioni illustrate nel paragrafo 5 2 5 e attribuendo al nodo polo peso doppio rispetto ai nodi satellite La soluzione del sistema di equazioni lineari ottenuto dall assemblaggio dei bilanci delle forze ai vari nodi stata ottenuta con un codice realizzato all uopo in linguaggio C utilizzante la libreria PETSc3 per la soluzione di sistemi di equa zioni lineari 5 2 5 La definizione delle propriet materiali I valori di modulo di elasticit E sono stati derivati dalle informazioni del datset CT tramite due operazioni e la relazione di calibrazione densitometrica del dataset CT che correla i va lori di Hounsfield Unit HU con valori di densit minerale ash density individuata in fase di acquisizione della tomografia tramite lo European Spine Phantom 82 Pash 0 0007851 HU 0 1439277 5 1 cm e la relazione tra densit minerale e modulo di elasticit assunta dalla lette ratura 83 E 10 563 p22 GPa 5 2 ash Il
141. ione temporale Si ottiene in questo modo il quadro completo della struttura logica utilizzata strumento molto importante ai fini tanto della rappresentazione del modello di analisi quanto per l implementa zione di procedure per il calcolo per via numerica della soluzione Un efficace rappresentazione sintetica del modello fisico approntato per la solu zione di un problema offerta dai diagrammi di Tonti si tratta di grafici che evi denziando il percorso che conduce alle equazioni risolutive passando attraverso le varie equazioni di struttura e costitutive Inoltre evidenziano la suddivisione in variabili di sorgente e configurazione e presentano le attribuzioni delle varia bili stesse agli enti geometrici e temporali Uno schema generale di riferimento presentato in figura 2 6 Il diagramma tridimensionale e si compone di due piani il secondo piano rappresentato a tratto fine ospita le grandezze integrali nel tempo in primo piano in grassetto sono invece collocate le grandezze che si ottenute dalle precedenti tramite un operazione di differenza nel tempo ossia l analogo discreto della derivazione temporale In ciascun piano figurano due colonne di caselle cui corrispondono i simboli P L S e V che si riferiscono rispettivamente ai punti alle linee alle superfici e ai volumi La colonna di sinistra rappresenta le entit primali quella di destra le corrispondenti entit duali Da notare che gli elementi duali in ordine ovvia m
142. ioni della costellazione locale Questi stessi valori riportati in scala logaritmica nel grafico di figura 4 17 evidenziano come l ordine di convergenza per le varie interpolazioni sia pari a 2 105 4 L approccio meshless con il Metodo delle Celle e lt 0 log 8 4 e 5 satelliti 8 satelliti lt lt T 1 4 1 2 1 0 8 0 6 0 4 0 2 log lato Figura 4 17 Errore medio rispetto al passo di griglia per le implementa zioni con 4 5 e 8 nodi satellite relativamente alla griglia regolare di punti va rilevata la sostanziale equivalenza delle soluzioni ottenute utilizzando 4 e 5 satelliti per le quali si ottengono gli stessi valori di scarto quadratico medio l uso di una costellazione composta di 8 satelliti registra invece a parit di ordine di convergenza ottenuto un aumento di precisione 4 3 2 Sistemi piani deformabili Nei paragrafi seguenti verr illustrata l estensione all ambito dello studio dei cor pi deformabili in regime di piccoli spostamenti dell approccio meshless ad inter polazione locale polinomiale descritto in precedenza Si consideri dunque un sistema elastico piano di geometria nota sia il materia le omogeneo e isotropo e caratterizzato da comportamento elastico lineare con modulo di elasticit E e modulo di contrazione trasversale v Il corpo sia assog gettato a un sistema di forze e di vincoli noti tali da generare uno stato piano di sollecitazione di sforzo oppure di de
143. ioni in 5 distinte condizioni di sollecitazione il test conserva ancora un valore puramente 149 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle indicativo avendo due sostanziali limitazioni e stato condotto su un singolo esemplare di femore normale sia sotto il profilo anatomico che sotto il profilo delle condizioni di mineralizzazione e stato realizzato costruendo il modello MCM su un dataset CT sottocam pionato perdendo con questo sia dettaglio di definizione geometrica del segmento osseo che accuratezza nella mappatura delle propriet materiali Affinch una procedura di modellazione numerica possa ambire a costituire ele mento orientativo nella pratica clinica imprescindibile la valutazione dell atten dibilit della precisione numerica ottenibile con un processo che prende il nome di validazione Due elementi focali di questo processo 85 sono la variabilit in tersoggettiva ossia la verifica del modello numerico quando applicato in casi distinti e il confronto con esperimenti controllati in vitro In questo capitolo verr dunque descritto uno studio di validazione della pro cedura di modellazione MCM La base del confronto stata una campagna di misure sperimentali precedentemente condotta su 8 distinti esemplari di femore 6 diverse configurazioni di sollecitazione in campo elastico sono state applicate su ciascun femore misurando le corrispondenti deforma
144. izionamento si trover comunque ad essere nodo di superficie in secondo luogo le misure estratte dal modello FEM mediano i valori calcolati ai centroidi di tutti gli elementi esterni ossia aventi una faccia sulla superficie afferenti al nodo questa operazione muove dalla considerazio ne che nella pratica sperimentale l estensimetro ha una griglia di rilevazione che 160 6 3 Validazione con le misure sperimentali copre una superficie circolare di diametro pari a circa 3 mm misura che inci dentalmente corrisponde alla dimensione media degli elementi superficiali del modello FEM Queste considerazioni hanno pertanto ispirato una strategia alternativa per il cal colo dei valori di deformazione in corrispondenza della posizione stimata de gli estensimetri Anzich considerare i valori principali di deformazioni calco lati nel punto geometrico corrispondente agli estensimetri si adotta la seguente procedura e si individua il voxel v in cui cade la posizione dell estensimetro e siindividua l insieme V dei voxel limitrofi il cui centroide disti dal centroide di v meno di un valore R ossia si individuano i voxel il cui baricentro cada entro una sfera di raggio R e di centro il baricentro di v e tra i voxel dell insieme Y si individua l insieme Y pg dei voxel caratterizzati dal fatto di avere valore medio di HU ai vertici maggiore della soglia HU I valori di deformazione principale associati alla posizione nominale
145. izione di un modello e la superfici tassellizzate tridimensionali un terzo approccio permette di elaborare direttamente in maniera automatica l intera pila di immagini con siderata quale dataset tridimensionale Questo metodo il cui antenato il Marching Cube 6 consente l estrazione automatica di superfici di iso densit dal dataset CT la superficie rappresentata da simplessi poligona li tipicamente triangoli che costituiscono una superficie denominata tas sellizzata Il difetto di questo metodo che produce spesso superfici con contraddizioni topologiche quali ad esempio triangoli mancanti o sovrap posti inoltre si viene a generare solitamente un numero molto elevato di simplessi che nel caso di un femore umano pu tranquillamente arrivare all ordine delle migliaia di elementi e le superfici matematiche tridimensionali l ultimo gruppo di tecniche di segmentazione fa uso di superfici matematiche elasticamente deformabi li la corrispondenza con il contorno dell organo ottenuta attraverso la minimizzazione di un opportuna funzione di costo definita sul campo sca lare del dataset CT Di fatto questo approccio ancora in fase di sviluppo e a tutt oggi non pienamente utilizzabile inoltre le implementazioni at tuali richiedono la definizione empirica di una serie di parametri per ogni singolo dataset e si pu pertanto affermare che l algoritmo non sia ancora automatico Agli effetti pratici attualm
146. l approccio voxel mesh Si tratta di una procedura completamente automatizzata che permette diretta mente la creazione di una mesh cartesiana appoggiata sulla struttura del data set CT La disposizione regolare di punti in cui sono noti i valori di assorbimen to della radiazione offre una naturale discretizzazione basata sui parallelepipedi costituiti dai voxel Esistono dei sofisticati algoritmi basati sul valore locale del gradiente del campo scalare di assorbimento della radiazione e dunque come si visto di densit per individuare gli elementi superficiali o per realizzare raffittimenti localizzati nelle zone di transizione fra tessuti differenti 1 3 2 2 L individuazione della superficie di contorno la segmentazione Un metodo alternativo comporta invece l individuazione di una superficie di con torno l operazione con cui in un immagine diagnostica si individua la superficie di contorno di un organo in genere denominata segmentazione Per realizzare questa operazione esistono dei codici commerciali che fanno uso di vari algorit mi pi o meno automatizzati In termini generali le tecniche di segmentazione possono essere classificate in quattro categorie fondamentali 15 1 La creazione di modelli numerici subject specific di segmenti ossei e le nuvole di punti il contorno definito da un insieme di punti nello spa zio questa categoria di fatto poco usata per le poche garanzie di continui t della superficie di co
147. l dataset con valore di HU gt 250 corrispondente ad un valore di densit minerale indicativamente assimila bile a quello dell osso corticale presente nella regione distale del femore e sulla superficie della regione prossimale E HU di figura 6 5 con lo spettro dei valori HU delle acquisizioni CT riporta to in figura 6 6 si noti come le relazioni utilizzate prevedano che E si annulli per valori negativi di HU al contempo l istogramma testimonia un elevatissimo numero di punti con HU prossimo allo zero Al valore HU 0 che corrisponde all acqua nel caso specifico propriamente il materiale in cui immerso il femore corrisponde un valore di E pari a circa 310 GPa Questo fatto ha dunque ori ginato una sezione resistente ampia attorno ad ogni femore contribuendo alla parziale distribuzione del carico e riducendo quindi le deformazioni puntuali ed aggregate nel segmento osseo Questa considerazione ha suggerito di modificare l equazione 6 4 facendo la valere solo al di sopra di un valore di soglia HU al di sotto del quale viene assunto per E un valore convenzionale trascurabile 6 5 TE a h a HU b k f se HU gt HU oyu se HU lt HU non nullo per evitare problemi di singolarit della matrice di rigidezza del problema 157 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle E GPa 25000 20000 15000 10000 5000 HU 0 J 0 500 1000 1500
148. lancio corretta priva di ambiguit e di errori 3 2 6 Testnumerico La precisione della soluzione numerica ottenibile con l implementazione CM pro posta stata valutata risolvendo un problema di test per il quale disponibile in letteratura una soluzione in forma chiusa Il problema trattato quello della sfera cava pressurizzata di grosso spessore pa ragonabile al diametro interno della cavit di materiale omogeneo e isotropo e a comportamento elastico lineare in regime di piccoli spostamenti La figura 3 7 riporta lo schema di riferimento in queste condizioni il problema ha simmetria sferica e il campo di spostamenti risulta radiale rispetto al centro della sfera ed definito dalla relazione Pa b u r E dd 1 2v r 1 v 572 3 68 essendo a e b i diametri interno ed esterno della sfera rispettivamente P il va lore della pressione interna E e v il modulo elastico e il modulo di Poisson del materiale r la distanza dal centro della sfera Il problema replicato completamente descritto dai valori riportati nella ta bella di figura 3 7 31 cap 7 68 3 2 Problemi di elasticita lineare tridimensionali parametro valore unit di misura a 1 m b 2 m E 200 GPa v 0 26 m P 10 m Figura 3 7 Sfera cava pressurizzata schema del problema e notazioni di riferimento Beneficiando della simmetria sferica del problema si considerata solamente una porzione della geometria corri
149. lazione all estensione degli enti geometrici di riferimento e distingue le variabili globali dalle variabili locali e le variabili globali chiamate anche variabili integrali sono funzioni di do minio e sono grandezze associate ad elementi discreti dello spazio e del tempo punti P linee L superfici S volumi V istanti I ed intervalli di tempo T Le variabili che vengono rilevate sperimentalmente sono di fatto delle va riabili globali e questo fatto legato alla metodologia stessa della misura che associa di necessit le grandezze fisiche ad elementi discreti dello spa zio cos che un contenuto di massa resta naturalmente associato ad un volume un flusso di calore ad una superficie un impulso ad un intervallo di tempo e le variabili locali note come funzioni di campo sono invece delle fun zioni di punto analiticamente dipendenti dal posto e dal tempo e sono le grandezze comunemente utilizzate dalla formulazione differenziale Esse sono ottenute dalle variabili globali a mezzo degli operatori di densit e di tasso ossia dal limite del rapporto fra il valore della grandezza globale e l estensione dell elemento spazio temporale cui la grandezza riferita al tendere a zero di questa estensione esse sono pertanto slegate dagli attributi geometrici di partenza 27 2 Il Metodo delle Celle Va infine precisato che le variabili globali associate ai punti quali ad esempio lo spostamento nella meccanic
150. li di deformazione calcolati ai cen trodi degli elementi del dal modello FEM replicante la configurazione di carico neutral del femore identificato come 2501 Essendo i due modelli registrati spazialmente rispetto al medesimo riferimen to la procedura di confronto ha previsto e individuazione delle coordinate del centroide dell elemento 167 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle e l individuazione del voxel del modello MCM corrispondente al centroide considerato e definizione della funzione di interpolazione trilineare dei valori delle com ponenti di spostamento ai vertici del voxel e calcolo dei valori principali del tensore di deformazione per derivazione della funzione di interpolazione Questa procedura non stata estesa a tutti gli elementi del modello FEM essen do stati dal confronto esclusi alcuni elementi per i quali l indagine non risultava essere significativa sono stati esclusi nel dettaglio e gli elementi con centroide distante meno di un diametro femorale dalla base di incastro diametro valutato all imposta dell incastro in questa zona si hanno infatti fenomeni locali non significativi ai fini del confronto verificata la coerenza globale tra i due modelli in termini di campo di spostamenti e gli elementi di superficie del canale endostale del femore in questa zona gli elementi del modello FEM hanno inferiore qualit posizionamento spesso
151. linea Orientazione esterna di un punto l orientazione interna del volume che contiene il punto a U 2 Figura 2 3 Orientazioni interna ed esterna degli elementi geometrici di base nello spazio tridimensionale della vite Osservazione Dall esposizione dei vari casi di orientazione esterna emerge una significativa dualita tra oggetti muniti di orientazione interna e oggetti dotati di un orienta zione esterna Questo legame si esprime nel fatto che un orientazione interna esterna di un oggetto induce un orientazione esterna interna sul suo duale Questa regola della mutua induzione ha una piccola eccezione i punti sono tra dizionalmente orientati internamente come pozzi questo indurrebbe di necessit un orientazione esterna dei volumi entit duali dei punti nello spazio tridimen sionale le cui facce verrebbero a ricevere tutte un verso di attraversamento po sitivo dall esterno all interno del volume Tradizionalmente per si considerano 32 2 3 L orientazione degli elementi temporali quasi sempre positive le normali alle facce uscenti dal volume assumendo quindi unorientazione discorde da quella suggerita dall orientazione interna di punti 2 3 L orientazione degli elementi temporali Alcune grandezze fisiche sono riferite al tempo Per le variabili locali il tempo una delle coordinate da cui la grandezza pu dipendere al pari delle coordi nate spaziali L espression
152. lizzazione i dati delle immagini CT vengono normalizzati e approssimati a valori interi Con questa operazione i valori floating point x y corrispondenti al valore del 13 1 La creazione di modelli numerici subject specific di segmenti ossei coefficiente di assorbimento lineare nel pixel di coordinata x y sono trasformati in valori interi utilizzando un operazione di normalizzazione quale x y Bo 1 1 CT x y 1000 HH20 dove umo il valore del coefficiente di attenuazione dell acqua I valori CT x y sono noti come Hounsfield Unit HU e rappresentano i valori in uscita di tutto il processo di elaborazione dei segnali registrati 1 3 2 La costruzione della mesh dai dati CT La discretizzazione della struttura ossea individuata dai dati CT con la creazio ne di una mesh tridimensionale una fase molto molto onerosa in termini di tempo nella definizione di un modello ad elementi finiti peraltro anche una fase molto delicata poich la qualit della mesh influenza la qualit dei risultati e per poter ottenere una buona accuratezza dei risultati la mesh deve soddisfare condizioni molto restrittive L adozione di un sistema automatico per la generazione della mesh una scelta obbligata per rendere un processo di analisi compatibile con i tempi della prati ca clinica ovvero quando si voglia condurre un analisi su un elevato numero di campioni del tutto evidente che una procedura manuale comportere
153. ll equazione di Laplace che relativamente al problema di elasticit lineare in 174 regime di piccoli spostamenti problema della trave di Timoshenko I tre me todi sono stati poi valutati nella prospettiva applicativa della realizzazione di un modello numerico automaticamente costruito sui dati CT l approccio basato sull interpolazione con funzioni di base radiale ha incontrato ostacoli di caratte re computazionale dovuti alla significativa complessit algebrica dell espressio ne analitica dell equazione di bilancio questo approccio stato quindi escluso non essendo compatibile con le risorse algoritmiche dei software di Computer Algebra System disponibili Parimenti escluso stato l approccio basato sull ap prossimazione mediante polinomi di grado parziale uniforme rispetto alle coor dinate il fatto di avere i nodi disposti su una griglia regolare produceva infatti malcondizionamenti nella matrice di interpolazione per bassi gradi del polino mio laddove per gradi pi elevati si realizzava un estrema complessit formale unitamente all asimmetria della distribuzione dei punti attorno ad ogni nodo di volta in volta esaminato La formulazione meshless basata sulla definizione di un complesso locale di celle primali si invece rilevata adatta e di implementazione immediata una mesh locale di tetraedri viene definita per ogni punto del dataset CT con i nodi primi vicini nell allineamento di griglia l attribuzione delle pr
154. lla seguente tabella 2 5 Un complesso di celle inteso come insieme piccole porzioni adiacenti che compongono un dominio formato da punti spigoli facce e volumi topolo gicamente questi elementi sono considerati celle di dimensione rispettivamente zero uno due e tre e si parla infatti di 0 celle 1 celle 2 celle e 3 celle Dal punto di vista topologico un complesso di celle non gi un insieme di volumi quanto piuttosto un insieme di celle di varia dimensione Dal punto di vista formale ad 36 2 5 I complessi di celle elementi del complesso primale elementi del complesso primale punti P oS volumi V linee L superfici S superfici S lt gt linee L volumi V lt gt punti P duale di ha duale di hi j duale di hi Figura 2 5 Celle primali e celle duali nello spazio a 2 e a 3 dimensioni una cella primale di ordine p corrisponde un elemento del sistema duale di ordi ne n p e viceversa essendo n la dimensione dello spazio in cui si considera il sistema di celle Altro aspetto di dualit riguarda l orientazione degli elementi dei due sistemi di celle la figura 2 3 evidenzia come un orientazione interna delle celle primali induca un orientazione esterna delle corrispondenti celle duali La figura 2 5 illustra la formazione di celle duali nel piano e nello spazio con diverse forme di celle primali Va rilevato il fatto che tanto la scelta del punto all interno di ogni cella primal
155. lle colonne associate ai nodi satellite nella numerazione globale dei nodi 8se M l indice associato ad un nodo le righe e le colonne associate hanno indice 2M 1 e 2M 91 4 L approccio meshless con il Metodo delle Celle 4 2 2 1 Testsulla trave di Timoshenko L indagine sulla precisione dell implementazione stata condotta relativamente al problema della trave di Timoshenko per il quale si dispone di una soluzione analitica di confronto Si tratta di una trave incastrata ad una estremit e as soggetta ad una distribuzione parabolica di sforzi di taglio all estremit opposta come illustrato in figura 4 10 y L i Figura 4 10 Trave di Timoshenko schema schema statico condizioni di carico riferimento cartesiano e notazioni adottate si noti come la condizione di vincolo sia approssimata e valevole solo ai fini delle stime e dei calcoli tecnico ingegneri stici prevedendo la soluzione teorica per x 0 spostamenti non nulli Per il problema relativamente ad un materiale omogeneo e isotropo il campo di spostamenti u v definibile in forma chiusa dalle relazioni i a ox 3x x 2 y r 4 13 2 PA pr da 65 92 essendo e F ev il modulo di elasticit e il modulo di Poisson rispettivamente del materiale e P la risultante degli sforzi di taglio e Le Dlalunghezza e l altezza della trave e zey le coordinate del punto nel riferimento di figura 4 10
156. lleli ai lati nel dominio sono stati posizionati dei punti secondo una griglia regolare parallela ai lati da questa griglia con spostamenti random e facendo in modo da avere un punto in ogni vertice e un certo numero di punti posizionati lungo ciascun lato si sono poi ottenute delle distribuzioni casuali di punti nel dominio Sui punti del contorno sono stati imposti i valori assunti dalla funzione armonica w x y e cosy 4 41 Come parametro di valutazione della precisione della soluzione numerica viene utilizzato il valore dello scarto quadratico medio rispetto alla soluzione teorica Ni ari ice Dili a wi 4 42 essendo N il numero di nodi interni quelli su cui non sono stati imposti i valo ri della funzione soluzione del problema w il valore assunto dalla soluzione teorica nel punto i e w il valore ottenuto dalla soluzione numerica nello stesso punto Iper i punti posizionati sui lati lo spostamento avviene lungo il lato stesso mentre i punti posizionati sui vertici non subiscono alterazioni di posizione 104 4 3 Interpolazione locale con funzioni polinomiali Il problema cos posto stato risolto con la procedura meshless ad interpola zione polinomiale sopra descritta implementata con riferimento a costellazio ni locali composte da diverso numero di nodi ossia con differenti polinomi di interpolazione locale Ty a bx cy dx ey Ts a bz cy dr ez
157. lli numerici subject specific di segmenti ossei densit degli otto punti della griglia CT che si trovano attorno al centroide del l elemento 19 Questi metodi sono di implementazione assai semplice ma pos sono produrre dei risultati inaccurati quando la dimensione degli elementi sia significativamente maggiore della spaziatura della grigia di punti CT Un secondo approccio proposto inizialmente con limitazioni 20 e succes sivamente reso disponibile in una formulazione generalizzata nel codice open source BONEMAT _ VI 21 assegna agli elementi la media dei valori associati ai punti del dataset CT che cadono all interno di ogni elemento In questo caso anche se la dimensione dell elemento significativamente maggiore della spazia tura della griglia CT il risultato ottenuto ragionevolmente accurato Tuttavia neppure questo metodo pu produrre risultati soddisfacenti quando la dimensio ne degli elementi sia al contrario pari o minore delle dimensioni dei voxel Per questa ragione sarebbe importante poter disporre di un codice in grado di trat tare entrambe le situazioni e che completamente automatizzato potesse essere utilizzato anche con mesh non strutturate Un codice con queste caratteristiche chiamato BONEMAT_V2 stato in effetti realizzato dai Laboratori di Tecnologia Medica dell Istituto Ortopedico Rizzoli di Bologna ed stato verificato per con fronto con il BONEMAT_V1 in tre differenti tipi di osso 1 3 3 1
158. lo di elasticit E con la particolare procedura di modellazione in cui per evitare ogni procedura di discriminazione tra tessuto osseo e ambiente circostante tutto il datset stato oggetto della procedura di cal colo ivi compresi i nodi collocati nell acqua di contenimento del femore durante l acquisizione CT Quello che di fatto venuto a crearsi che in virt dell adottato legame HU E i nodi collocati nel materiale circostante il femore hanno ereditato una resistenza meccanica che pure di modesta entit ha contribuito a creare un sostegno resi stente circostante in grado di riprendere parte del carico applicato Il fenomeno suggerito dall immagine di figura 6 4 che rappresenta il campo di HU di uno dei femori e in cui sono stati selezionati con colorazione opaca i punti con HU superiore a 250 valore corrispondente a valori di densit minerale caratteristi ci dell osso corticale l immagine d quindi una rappresentazione approssimata della geometria ossea indicando come per essa sia circondata da altro materiale nella particolare descrizione materiale adottata questo materiale evidentemente abbondante assume una resistenza meccanica nel complesso non trascurabile Questo fatto stato reso evidente dalla comparazione della funzione aggregata 156 6 2 Le criticita emerse a Figura 6 4 Rappresentazione del campo HU del datset CT del femore denomi nato 2501 in opaco sono raffigurate le regioni de
159. lo nella numerazione globa le Le colonne in cui gli elementi di kp vanno collocati corrispondono alle colonne di indice pari al numero identificativo di ogni singolo nodo della costellazione nella numerazione globale 4 4 1 3 Alcuni test numerici Al fine di verificarne l accuratezza l approccio descritto stato applicato alla so luzione della equazione di Laplace V w 0 su un domino piano quadrato di lato unitario centrato nell origine di un sistema di riferimento cartesiano con assi paralleli ai lati nel dominio sono stati posizionati dei punti secondo una griglia regolare facendo in modo di avere un punto in ogni vertice e un certo numero di punti posizionati lungo ciascun lato Sui punti del contorno sono stati imposti i valori assunti dalla funzione armonica wy du e Cosy 4 83 120 4 4 Interpolazione locale con funzioni di base radiale soluzione del problema Come parametro di valutazione della precisione della soluzione numerica viene utilizzato il valore N ae de Da ct wi 4 84 essendo N il numero di nodi interni quelli su cui non sono stati imposti i valo ri della funzione soluzione del problema w il valore assunto dalla soluzione teorica nel punto i e w il valore ottenuto dalla soluzione numerica nello stesso punto Il problema cos posto stato risolto con la procedura meshless con interpolazio ne RBF sopra descritta implementata con riferimento a quattro diverse funzioni
160. losest point 88 risultato variabile per i vari modelli tra 0 5 e 0 9 mm essendo l osso in prima approssimazione un materiale poroso la densit apparente esprime il rapporto massa su volume totale di un campione valutata per ignizione di un campione valutando il peso della frazione residua 154 6 1 L esperimento di riferimento 6 1 4 Imodelli MCM I modelli MCM sono stati costruiti direttamente a partire dai datset delle acquisi zioni CT di ogni singolo femore Registrato il riferimento numerico con quello sperimentale sono state replicate le condizioni di carico e di vincolo di ogni singola prova il dataset a risoluzione originaria stato croppato a livello del cemento eliminando le informazioni a quota inferiore a quelle dell imposta della base di PMMA e a tutti i punti della superficie di base sono stati imposti spostamenti nulli replicando il vincolo spe rimentale di incastro Le forze sono state applicate con le componenti ottenute dalla registrazione spaziale al nodo del dataset pi prossimo al punto di appli cazione del carico sperimentale L algoritmo MCM stato applicato costruendo attorno ad ogni nodo una mesh locale di tetraedri a 4 nodi appoggiati ai nodi prossimi secondo le direzioni di allineamento della griglia del dataset CT Ai nodi sono stati calcolati i valori di E ottenuti a partire dai corrispondenti valori di HU con le relazioni riportate nel paragrafo 6 1 3 La dimensione dei
161. me pena l errato com puto delle forze agenti Per un nodo interno la regione tributaria si estender appieno attorno al nodo se il nodo giace invece sul contorno del solido in esame la regione tributaria sar 126 4 4 Interpolazione locale con funzioni di base radiale la sola parte del cerchio contenuta nel solido A y Figura 4 24 Costellazione e regione tributaria per il bilancio delle forze al generico nodo polo Pp Il calcolo delle forze di superficie agenti su ogni lato della regione tributaria ri sulta cos semplificato La forza F agente su una superficie orientata A espressa dalla relazione di Cauchy F f DndA 4 97 A essendo n il versore normale alla superficie e X il tensore degli sforzi qui inteso nella notazione matriciale da i al 4 98 Tay Oy Con riferimento alla figura 4 24 si ha infatti essendo t lo spessore del sistema piano che per il lato HK che il vettore area corrispondente risulta essere S 2td 0 da cui d Fruk ox d y dy d 4 99 d Fy aK Try d y dy d 127 4 L approccio meshless con il Metodo delle Celle da cui si ottiene 8 3d x d y 3d x3 9d x 10d 3 d x 4d 3d A 3 iv 8 d x 2d xi 2d d Vi dy B wi Fr HK TA F ug 8d 1 x 2dz 2d Yi yi C Wi 8 3d x d y 3d x3 9d x 10d 3 d x 32d 3d C ee 3 Wiv 4 100
162. me flusso della quantit di moto e per tanto l equazione di equilibrio pu essere vista come equazione di stazionariet del flusso della quantit di moto Appare infine rilevante accennare al fatto che il bilancio pi che una legge fisica in s un costrutto della pensiero logico dell uomo che trova applicazioni in molti ambiti applicativi La sua importanza oltre che concettuale quella di porre in relazione quantit accumulate prodotte transitate definite da altre leggi fisiche 2 6 2 Equazioni costitutive Le equazioni costitutive sono delle relazioni che definiscono la costituzione del sistema riportando i parametri caratteristici del mezzo in cui hanno luogo i fe nomeni studiati Queste equazioni caratterizzano il comportamento dei materiali presenti rispetto al problema studiato e pertanto identificano un modello di com portamento Sono esempi di legge costitutiva le leggi nella meccanica dei solidi che descrivono il modello di comportamento elastico plastico viscoelastico di un materiale Le equazioni costitutive sono anche chiamate equazioni di comportamento o nella termostatica equazioni di stato Caratteristiche essenziali delle equazioni costitu tive sono e legare variabili differenti variabili di configurazione con variabili di sor gente e contenere costanti materiali e essere valide in regioni a campo uniforme intendendo per campo uniforme l uniformit del gradiente del potenziale del campo
163. metodo pi accurato usando tutte le informazioni disponibili per inter polare i valori di HU all interno del dominio CT In questo caso la dimensione dell elemento finito non influenza l accuratezza nella stima di HU dato che la descrizione del campo scalare delle HU non uniforme anche all interno di ogni singolo voxel e tiene conto della distribuzione spaziale delle HU negli otto nodi circostanti i vertici della griglia di punti CT Calibrazione del data set CT I valori numerici prodotti da un acquisizione CT dipendono da molti fattori legati allo specifico esame tra cui i parametri fisici impostati dal radiologo come il KVP kVolt Peak o l X ray Tube Current la calibrazione realizzata in genere con l utilizzo di una campione di controllo o fantoccio necessaria per trasformare correttamente i valori misurati in valori significativi di densit di tessuto osseo Si assume sia valida una relazione di tipo lineare e che ogni utente possa adottare la propria relazione con i parametri opportuni L equazione di calibrazione allora Pa a BHU 1 3 dove p la densit assegnata in modo uniforme all elemento n della mesh HU il valore uniforme del numero CT a e 8 sono i coefficienti di calibrazione forniti dall utente Definizione del modulo di Young In letteratura si possono rintracciare molte relazioni che legano la densit del tessuto osseo al modulo di Young La relazione viene indicata dall utente e pu essere s
164. minio in esame sola mente a livello di punti nodi posizionati dentro la geometria e sul suo contor no come suggerito dalla figura 4 1 L approssimazione di una funzione scalare avviene tramite una relazione del tipo u x t 9 x ur t 4 1 Tes essendo con riferimento alle notazioni di figura 4 1 e uzil valore al nodo I di posizione x e Q R la funzione di forma associate al nodo J solitamente si tratta di una funzione non nulle solo entro una regione di influenza del nodo 2 definita da un valore d detto supporto del nodo J e S l insieme dei nodi per i quali la funzione di forma x non nulla L espressione formale dunque analoga a quella nota per gli Elementi Finiti con la differenza che ora le funzioni dell equazione 4 1 sono approssimanti e non gi interpolanti verificandosiu xr Questo fatto ha come di seguito specificato un importante ricaduta nell imposizione delle condizioni al contorno di tipo es senziale I vari metodi differiscono tra loro essenzialmente per e il tipo di approssimazione che viene effettuata si distinguono a tale pro posito i metodi intrinsechi SPH MLS dai metodi detti estrinsechi PUFEM hp clouds e il tipo di integrazione della forma differenziale in particolare alcuni metodi utilizzano la forma forte avvalendosi del metodo di collocazione per la costruzione del sistema lineare risolvente altri usano invece il metodo di Galerkin integran
165. modelli risultata variabile tra 2998800 e 3764640 nodi Il solutore stato costruito con un programma ad hoc in linguaggio C con chia mata a funzioni della libreria PETSc per l esecuzione su architetture parallele la soluzione del sistema lineare risolvente avvenuta con il metodo del Gradiente Coniugato con tolleranza di ciclo pari a le 7 6 1 5 Definizione dell accuratezza L elaborazione delle misure estensimetriche restituisce per ogni estensimetro e per ogni condizione di carico i valori principali di deformazione sulla superficie dell osso Questi valori sono stati usati come elemento di confronto per i valori calcolati con il modello numerico nei punti corrispondenti alla posizione degli estensimetri L assunzione di stato piano di stress sulla superficie alla base della procedu ra con cui sono stati determinati i valori principali di deformazione nel modello FEM come descritto nel paragrafo 5 2 6 sono assunti come valori di confronto i valori principali di deformazioni associati alle direzioni principali di tensione cui corrispondono i due valori di stress in valore assoluto pi distanti dallo 0 Per il modello MCM invece non vi sono assunzioni sulla definizione di superfi cie e il segmento osseo viene considerato un continuum con il suo intorno ma nipolando solamente valori di densit e di modulo di elasticit mappati sulla griglia regolare di punti prodotta dall acquisizione CT Il calcolo del tensore di d
166. modita implementativa verra adottato un sistema di riferimento locale ot tenuto traslando il sistema globale in modo tale da posizionarne l origine nel 116 4 4 Interpolazione locale con funzioni di base radiale nodo polo Come funzione di base radiale r si adotti a scopo puramente esemplificativo la funzione g r r 4 67 essendo r la distanza euclidea rispetto al nodo polo Le condizioni di interpolazione permetteranno di esprimere i pesi w in termini di temperature nodali t avendosi cfr equazione 4 64 il sistema di equazioni lineari Bewe Te 4 68 da cui we Te 4 69 La funzione approssimante appoggiata ai valori nodali incogniti delle tem peratura e i coefficienti wc dell interpolazione sono espressi a meno della matrice c che contiene le informazioni sulla disposizione geometrica dei punti della co stellazione in termini del vettore delle temperature nodali Te La procedura come per il caso dell interpolazione polinomiale esprimera il vet tore gradiente derivando la funzione T P ed esprimer il vettore flusso di calore unitario tramite la legge costitutiva del materiale infine definita una regione tri butaria attorno al nodo polo il flusso di calore uscente attraverso la sua superficie di contorno sar calcolato con un integrale di superficie della funzione che espri me il vettore flusso unitario Questo processo calcola un valore scalare il flusso netto di cal
167. modulo di Poisson stato assunto pari a 0 3 84 5 2 6 Determinazione dell accuratezza dei modelli numerici Le misure risultanti 26 sono state utilizzate per la validazione delle deformazio ni principali calcolate dai modelli numerici Shttp www mcs anl gov petsc petsc as 144 5 2 Il test esplorativo 5 2 6 1 Modello FEM Per il modello FEM il calcolo dei valori principali di deformazione avvenuto individuando le due direzioni principali di tensione cui corrispondessero valori di tensione pi lontani dallo zero ossia escludendo la direzione principale con minore valore di tensione associato Si suppone infatti uno stato di stress pia no sulla superficie condizione usualmente verificata nella simulazione numerica dove una delle tensioni principali calcolate solitamente assai prossima allo zero Si calcolano a questo punto le deformazioni principali corrispondenti da utiliz zare per il confronto con le misure sperimentali Va sottolineato come esse non corrispondano necessariamente ai valori principali di deformazione minimo e massimo poich contrariamente allo stato di tensione non possibile ipotizzare l esistenza di uno stato di deformazione piano in questa condizione i valori di deformazione principale minimo o massimo potrebbero corrispondere a direzio ni non tangenti alla superficie Nella fase di preparazione della mesh stata posta attenzione a che un nodo cor rispondesse alla posizione di
168. n eleva ta corrispondenza con i valori sperimentali verificata l indipendenza degli scarti 161 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle 3000 UE 2000 1000 5 olim Ateo I I Ai et M as ci exp DM FEM D MCM i R2 E R3 1000 2000 1000 5 EXP Fem O mcm MR2 E R3 1500 Figura 6 8 Valori di deformazione nel femore 2501 per la condizione di carico neutral per ogni deformazione sono riportate nell ordine il valore sperimenta le la predizione FEM la predizione MCM puntuale le predizioni MCM mediate con raggio pari rispettivamente a 3 e a 5 mm fra predizioni e misure sia dalle condizioni di carico che dall esemplare di femore ANOVA fattoriale p value gt 0 05 predizioni e misure sono state analizzate nel loro insieme valutando i parametri di accuratezza a livello aggregato La retta di regressione delle predizioni MCM rispetto alle misure sperimentali con pendenza prossima all unit intercetta prossima allo zero e coefficiente R pari a 0 94 manifesta l elevata correlazione tra predizioni e misure sperimentali Gli indicatori di accuratezza globale sono risultati essere soddisfacenti essendo 162 6 3 Validazione con le misure sperimentali il valore di errore massimo pari al 47 77 e l errore quadratico medio apri al 7 92 Come riportato in tabella 6 1 questi elementi di valutazione dell accu
169. n ampio range di valori di densit 8 9 Il passaggio fra le CT units e la densit di un tessuto biologico chiamato calibrazione del dataset Inoltre sono state individuate sperimentalmente delle buone relazioni empiriche tra densit e propriet meccaniche dei tessuti ossei 10 11 12 quindi in li nea di principio possibile derivare dai dati CT una distribuzione non omogenea delle propriet dei tessuti e tenerne quindi conto nel modello ad elementi fini ti subject specific In letteratura sono proposti vari metodi per realizzare questa operazione caratterizzati da differenti livelli di automaticit e dipendenti dalla tecnica che si adottata per creare la mesh Nel caso cos detto voxel mesh 13 14 il compito semplice il data set CT rap presenta un volume campionato nei punti di una griglia regolare e gli elementi della mesh cartesiana generata automaticamente con un criterio di sogliatura au tomatica hanno forma di parallelepipedo e sono ottenuti dall unione di un pre fissato numero di voxels Il trasferimento dei dati CT agli elementi della mesh dunque naturale potendo si associare direttamente ad ogni elemento il valore medio di densit radiologica dei voxel che lo compongono il passaggio al modulo di Young poi possibile utilizzando un opportuna relazione empirica riportata in letteratura Va sotto lineato che per la natura delle relazioni che legano valori CT alla densit e la 18 1 3 La de
170. na modello a elementi finiti subject specific 23 1 La creazione di modelli numerici subject specific di segmenti ossei 24 Capitolo 2 Il Metodo delle Celle Il Metodo delle Celle un metodo numerico basato sulla formulazione finita di retta delle equazioni dei campi fisici Attualmente i metodi numerici partono dalle equazioni di campo scritte in forma differenziale ottenute dall analisi di porzioni infinitesime dei sistemi fisi ci La scelta di fissare l attenzione su porzioni di spazio di estensione tendente a zero risponde all intento di ottenere una soluzione esatta dei problemi descritti questo obbliga a ridurre le grandezze fisiche che nascono associate ad intervalli di tempo o a elementi geometrici di estensione finita a funzioni del tempo e del punto in modo tale da poter applicare gli operatori matematici del calcolo diffe renziale Questa strategia costituisce certamente un importante struttura logica di studio ma comporta alcune difficolt in primo luogo le funzioni che defi niscono le grandezze fisiche devono essere derivabili rispetto alle coordinate di riferimento questo fatto comporta difficolt nella presenza di discontinuit geo metriche e nelle propriet dei materiali In secondo luogo la presenza di sorgenti concentrate astrazione che rappresenta sorgenti distribuite su estensioni piccole rispetto alle dimensioni del sistema introduce singolarit di carattere puramente matematico In ultima anali
171. ndini in regime di grandi deformazioni modello poroelastico delle cartilagini articolari evidenziando i potenziali vantaggi del l utilizzo del metodo meshless nelle simulazioni biomeccaniche Ciononostante manca ancora una applicazione a livello di organo con un modello tridimensio nale di un segmento osseo intero Lee et al 59 applicano una implementazione della Moving Least Square Ap proximation Technique per analizzare il comportamento meccanico di campioni cubici di osso trabecolare sottoposti a compressione con un modello tridimensio nale dei campioni sperimentali derivato da immagini u CT il lavoro dimostra le notevoli potenzialit rappresentate dall applicazione diretta del metodo meshless ai dataset di immagini cliniche per la costruzione di modelli numerici 4 2 Meshlessa celle locali con il CM Come stato evidenziato il Metodo delle Celle scrive delle relazioni algebriche tra variabili di sorgente o di configurazione che siano di natura globale Gli elementi geometrici e temporali di supporto per riferire queste grandez ze sono costituiti dalle celle di ordine n dei due complessi primale e duale Tra i due complessi di celle sussiste una relazione di dualit che si manifesta nella corrispondenza biunivoca tra gli elementi corrispondenti La chiave di volta dell approccio pu essere vista nella scrittura di un equazio ne di bilancio riferita alle celle duali Le celle duali sono ottenute da una suddivi sione
172. nfine fra la porzione di cella associata al nodo e il resto della cella Con la seconda strada si costruisce invece un sistema lineare locale che mette in relazione per ogni cella c esaminata le temperature nodali con i valori di calore generato assorbito dalla porzione di cella associata ad ogni suo nodo k T Q 3 33 questo sistema deve essere assemblato opportunamente sommando i termini nelle righe e per la matrice k anche nelle colonne del sistema risolvente globale di indice corrispondente ai numeri che identificano ogni nodo in una numerazio ne globale Un ultima considerazione riguarda i segni dei flussi di calore attraverso le su perfici essi sono strettamente correlati all orientazione esterna assunta per tali elementi di superficie Va qui ricordato fatto non solo formale bens sostanziale per l implementazione del Metodo che si deve provvedere ad orientare inter namente gli elementi del complesso di celle primale poich questa orientazione indurr automaticamente e in modo non equivoco l orientazione esterna degli elementi corrispondenti del complesso di celle duale 3 1 2 Test di convergenza Per realizzare un test di convergenza della soluzione numerica ottenibile con sim plessi nello spazio con interpolazione quadratica stata considerata l equazione di Laplace AT x y z 0 e due funzioni di riferimento adottate come test per la stima dell accuratezza della soluzione T x y z e siny e
173. nsion of the meshless approa ch for the cell method to three dimensional numerical integration of di screte conservation laws International Journal for Computational Methods in Engineering Science and Mechanics 7 2 69 79 11 2006 67 Charles A Micchelli Interpolation of scattered data Distance matrices and conditionally positive definite functions Constructive Approximation 2 1 11 22 December 1986 68 M J D Powell Radial basis functions for multivariable interpolation a review pages 143 167 1987 69 I S Raju D R Phillips and T Krishnamurthy A radial basis function approach in the meshless local petrov galerkin method for euler bernoulli beam problems Computational Mechanics 34 6 464 474 November 2004 70 X Zhang K Z Song M W Lu and X Liu Meshless methods based on collocation with radial basis functions Computational Mechanics 26 4 333 343 October 2000 71 Mehdi Tatari and Mehdi Dehghan A method for solving partial differen tial equations via radial basis functions Application to the heat equation Engineering Analysis with Boundary Elements 34 3 206 212 March 2010 72 Martin D Buhmann Radial Basis Functions Theory and Implementations Cambridge University Press 2003 ISBN 978 0 521 63338 3 73 M D Buhmann Radial basis functions Acta Numerica 9 1 1 38 2000 74 Fr d ric Magoul s Luis A Diago and Ichiro Hagiwara Efficient precon ditioning for image reconstruct
174. ntato in ternamente quando siano fissate delle orientazioni compatibili sulle sue facce di contorno 2 2 2 Orientazione esterna degli elementi spaziali L orientazione esterna di un elemento geometrico dipende dalla dimensione del lo spazio in cui l elemento considerato immerso Mentre l orientazione interna una propriet intrinseca l orientazione esterna cambia a seconda che l elemen to geometrico sia considerato nel piano o nello spazio Questo tipo di orienta zione meno immediata di quella interna e merita una pertanto una disamina dettagliata essendo il punto un entit geometrica priva di estensione lecito ritenere che una sua orienta zione sia priva di senso geometrico Anche se la nozione non ha un significato immediato essa ha comunque senso per le propriet formali che permette di mantenere Per analogia si pu pensare all elevamento a potenza scrivere n un espressione in s priva di significato se rapportata alla definizione che prevede di moltiplicare n per se stesso 0 volte La nota regola n n n estesa al caso in cui a b permette per di dare significato all espressione n che diviene come noto n n n n 1 Attribuire significato all espressione n permette di conservare le propriet formali dell elevamento a potenza 30 2 2 L orientazione degli elementi spaziali 2 2 2 1 Orientazione esterna di una linea La figura 2 2 rappresenta le possibili orientazion
175. nte il prodotto del tensore o per il versore normale alla superficie n Questo processo calcola due valori scalari le componenti della forza netta agente sulla superficie 00 p di confine della regione tributaria Qp del nodo polo P elabo rando le funzioni che esprimono i campi locali delle componenti di spostamento nell intorno del polo poich gli operatori matematici che intervengono in que sto processo derivazione operazioni algebriche e integrazione sono lineari ne consegue che l intero processo anch esso lineare Se dunque la funzione u x y viene espressa come somma di funzioni componenti u la forza totale Fan pu venire calcolata sommando i valori F o ottenuti applicando il processo ad ogni singola funzione addenda u Alla luce di questa considerazione tratteremo di seguito una sola delle funzioni che compongono la funzione approssimante i 0 Y Win 1 r2 wu 1 z zi y 7 ogg c 1 e xi y yi PE 2 vi ty Wine 1 1 wwe 1 z zi y yi centrata nel nodi generico i della costellazione C 4 4 2 2 Deformazioni Le componenti del tensore delle deformazioni vengono ricavate direttamente tramite le opportune derivate delle componenti dello spostamento Ou Ox le a yf dy ia 1 u dv 2 5 i x 4 94 4 x qi y yi a xi 1 Wiu _ Ausy yy a 1 win 4 z x yi Sr xi 1 Wiut 4 L approccio meshless con
176. ntorno teoricamente esisterebbero dei metodi tipici dei programmi CAD per approssimare la nuvola di punti con superfici ma tematicamente definite ma di fatto forse per causa della dispersione di pre cisione nel posizionamento dei punti non si a conoscenza di applicazioni in ambito biomeccanico le pile di contorni una seconda categoria di algoritmi elabora separata mente ogni singola immagine della pila prodotta dalla CT per ogni im magine viene utilizzato un algoritmo di segmentazione quale il conjugate gradient follower 3 e il contorno viene cos definito da una serie di curve piane Tutti i metodi convenzionali basati su derivative methods sono per al pi semi automatici anche nel caso della segmentazione di un os so facilitato da immagini fortemente contrastate un intervento manuale si rende in genere necessario nelle regioni epifisarie 4 Segmentate tutte le Dataset TC Figura 1 4 Segmentazione delle immagini e produzione dei contorni per ciascu na delle immagini del data set CT immagini l insieme dei singoli contorni pu essere interpolato per ottenere una superficie chiusa questa operazione sempre fattibile adottando come superficie un insieme di simplessi 5 ma se si vuole fare uso di superfici matematicamente descritte da polinomi di ordine superiore si deve in ge nere provvedere alla correzione manuale dei coefficienti del polinomio per ogni singolo contorno 16 1 3 La defin
177. num vit 4 14 N avendo indicato con e N il numero dei nodi interni dunque non vincolati e per cui ha senso il confronto con la soluzione teorica e con u e con v le componenti del vettore spostamento secondo la direzione x e y rispettivamente del sistema di riferimento globale e con il pedice num un valore di soluzione numerica e con il pedice un valore di soluzione teorica Il calcolo stato condotto distribuendo dei punti all interno del dominio secondo una griglia regolare equispaziata i valori di 6 al variare del passo di griglia sono riportati in tabella 4 3 e in scala logaritmica in figura 4 11 4 3 Interpolazione locale con funzioni polinomiali Nel precedente paragrafo stata presentato un approccio meshless immediata mente offerto dalla struttura concettuale e procedurale del Metodo delle Celle 93 4 L approccio meshless con il Metodo delle Celle parametro valore unit di misura lato L 8 m 1 0 0082415 D 2 m 0 5 0 0071274 E 1e3 GPa 0 25 0 0066124 v 0 25 0 125 0 0063476 P 2 N Tabella 4 2 Parametri descrit tivi il problema della trave di Ti moshenko usato per il confronto Tabella 4 3 Valore dellerrore medio 6 per i punti interni al dominio al variare del passo della numerico griglia di punti 2 08 5 e 2 1 2 124 wD 214 4 o hae e 2 16 4 2 18 5 e 2 2 T T T T I 1 0 8 0 6 0 4 0 2 0 log lato Figura 4 11 Test
178. nzione multiquadratica sono dette funzioni loca lizzate per esprimere il fatto che tendono a zero al tendere dell argomento a infinito il comportamento contrario r 0 per r oo della funzione multi quadratica viene descritto con l appellativo non locale 114 4 4 Interpolazione locale con funzioni di base radiale Assunta una funzione per determinare i valori degli N coefficienti w si utilizza la condizione di interpolazione Va Ei 4 63 ottenendo il sistema lineare r d r12 113 rin Wi di d 721 d 722 723 72 w2 dy P 731 O r3 2 r33 Esn 4 wae 4 ds 4 64 rma Plna Alna o Olan lun Xda avendo inteso indicare r x x con evidente significato della notazione il sistema pu scriversi in forma sintetica come w d 4 65 La matrice una matrice quadrata di dimensione NxN ed detta matrice di interpolazione Si pu dimostrare che esiste una una classe C di funzioni di base radiale che include la funzione gaussiana e le funzioni multiquadratiche per cui se i punti di interpolazione x1 2 p R sono distinti la matrice di interpolazione non singolare Inoltre per la funzione gaussiana e per la funzione multiquadratica inversa la matrice di interpolazione risulta anche essere definita positiva Dal punto di vista dell applicazione qui trattata i principali vantaggi dell ado zione delle RBF per interpolare lo
179. o Energy X ray Absorptiometry ha collocato questi campioni tra la condizione di osteopenia e quello di severe osteoporosis Precedentemente alle pro ve meccaniche non distruttive stata realizzata un acquisizione CT di ciascun esemplare secondo le specifiche di seguito riportate 6 1 1 Le misure sperimentali Ciascun femore stato cementato alla base nella regione condilare in un conteni tore di acciaio a sezione rettangolare per mezzo di PMMA Sulla superficie sono stati posizionati 15 estensimetri secondo ciascun aspetto anatomico anteriore laterale posteriore e mediale su 4 livelli collocati principalmente nella zona epi fisaria come illustrato in figura 6 1 La superficie di applicazione di ciascuno dei sensori stata preparata secondo una procedura validata in letteratura 78 per le prove su campioni idratati di osso cadaverico Le 15 rosette estensimetri che applicate erano del modello KFG 3 120 D1711L3M2S Kyowa Tokyo Japan aventi una lunghezza di griglia di 3 mm La tensione applicata stata limitata a 0 5 V per griglia per evitare effetti di surriscaldamento Le deformazioni sono state acquisite ad una frequenza di 10 Hz con un filtro passa basso in frequenza fissato a 1 Hz 130 mm Figura 6 1 Rappresentazione dell aspetto anatomico anteriore e laterale di un femore destro strumentato la posizione nominale degli estensimetri indicata relativamente alla lunghezza biomeccanica del femo
180. o B Calvo M A Martinez J M Garcia and J Cegonino On the employ of meshless methods in biomechanics Computer Methods in Applied Mechanics and Engineering 194 6 8 801 821 2005 59 James D Lee Youping Chen Xiaowei Zeng Azim Eskandarian and Morton Oskard Modeling and simulation of osteoporosis and fracture of trabecular bone by meshless method International Journal of Engineering Science 45 2 8 329 338 2007 60 Luigino Zovatto and Matteo Nicolini Improving the convergence order of the meshless approach for the cell method for numerical integration of di screte conservation laws International Journal for Computational Methods in Engineering Science and Mechanics 8 5 273 281 september 2007 61 Enzo Tonti and Federica Zarantonello Algebraic formulation of elastostatics the cell method CMES 39 3 201 236 2009 62 S P Timoshenko and N Godier Theory of Elasticity McGraw Hill 3rd edition 1987 63 Charles E Augarde and Andrew J Deeks The use of timoshenko s exact so lution for a cantilever beam in adaptive analysis Finite Elements in Analysis and Design 44 9 10 595 601 2008 64 O C Zienkiewicz and Robert L Taylor The Finite Element Method Basic Con cepts and Linear Applications volume 1 McGraw Hill London iv edition 1989 65 Valeriano Comincioli Analisi Numerica Metodi Modelli Applicazioni MCGRAW HILL 1990 183 66 Luigino Zovatto and Matteo Nicolini Exte
181. o MCM non discrimina in maniera dicotomica l osso dal non osso ma piut tosto mappa delle propriet materiali la densit il modulo di elasticit sui punti della griglia prodotta dal dataset CT vedendo lo spazio come un continuum descritto puntualmente e il confronto puntuale mette in relazione informazioni puntuali estratte da due entit geometriche il tetraedro e il voxel di dimensioni ed estensione spesso molto differenti tra loro se gli elementi del tetraedro in superficie hanno lato medio di 3 mm e risultano a forma regolare il voxel da cui vie ne estratta l informazione di deformazione nel modello MCM ha lato pari a 0 59 mm e altezza variabile tra 1 3 mm nella regione prossimale e 5 mm nella regione dello stelo Questo fatto rende peraltro impossibile attua re la strategia di mediazione sui voxel limitrofi descritta per il confronto con i dati sperimentali poich non sussiste in termini generali conformit geometrica tale da poter guidare la selezione dei voxel prossimi Queste considerazioni circa la natura locale delle differenze nelle predizioni del le deformazioni sono avvalorate dal confronto del modulo degli spostamenti ai nodi illustrato in termini relativi nel grafico di figura 6 14 nel quale si evince come i due modelli siano coerenti 170 6 4 Comparazione sistematica con modello FEM di riferimento 100 T 0 500 el 80 7 4000 60 7 3000 40 2000 20 4 1000
182. o al quale gli au tori illustrano la validazione dell analisi con misure sperimentali condotte su un cervello suino Sempre nell ambito dei modelli soft tissue Pena et al 56 applicano il NEM al problema delle grandi deformazioni in materiali non isotropi fibro rinforzati gli autori illustrano un modello di legamento umano confrontando i risultati con quelli ottenuti da una analisi FEM accoppiata Nello specifico ambito della biomeccanica dei tessuti ossei l uso di approcci meshless invece molto limitato solo tre lavori in letteratura riportano l utilizzo di algoritmi meshless nella predizione del comportamento meccanico di segmen ti ossei Liew et al 57 si avvalgono del Reproducing Kernel Particle Method per analizzare la distribuzione delle deformazioni nella parte prossimale del fe more umano in una simulazione dell effetto dell invecchiamento e della necrosi di tessuti oltre all interesse costituito dall essere un primo lavoro che utilizza una metodologia meshless nello specifico ambito biomeccanico l adozione di un mo dello semplificato bidimensionale e l adozione di un materiale base omogeneo costituiscono un significativo limite nella prospettiva delle applicazioni cliniche 78 4 2 Meshless a celle locali con il CM Lo studio di Doblar et al 58 esplora invece l applicazione del NEM in una relativamente ampia casistica rimodellamento adattivo del tessuto osseo com portamento iperelastico dei te
183. o cartesiano locale con origine in questo nodo e assi paralleli al sistema di riferimento globale Tra i nodi prossimi al nodo polo si selezioni una costel lazione di nodi satellite S per semplicit di illustrazione si assuma ora che la costellazione sia composta da 4 satelliti ofr 64 cap 7 96 4 3 Interpolazione locale con funzioni polinomiali Pe Ne Te ANTE OE LTR 5 y 2 S4 P i o R x l as E i Sa S3 g e e S2 Oe A E E gA IAA PI Figura 4 12 Costellazione di 4 satelliti attorno ad un generico nodo in terno assunto come polo Nell intorno del polo P il campo scalare studiato verr approssimato inter polando i valori T ai nodi della costellazione polo P pi satelliti S con una funzione polinomiale di grado due T P T y a br cyt dz ey 4 15 I quattro coefficienti sono determinabili imponendo che il polinomio assuma in ciascun punto della costellazione il valore scalare assunto dal campo in quel punto T P T ziyi a bzi cyi dr ey Ty 4 16 Le componenti del vettore gradiente della temperatura gc in un punto P x y sono esprimibili in forma analitica derivando la funzione di interpolazione locale _ ar Ju dr OT Fy c 2fy b 2dx 4 17 Il flusso unitario di calore q correlato con il vettore gradiente delle temperature a meno della conducibilit termica k del mezzo qe kgc 4 18
184. o dall estrema onerosita in termini di tempo e di risorse qualificate necessarie per la preparazione di ogni singolo modello la fase di pre processing scontornamento del dataset CT definizione matematica della superfici di contorno dell osso generazione della mesh mappatura delle propriet mate riali solo in parte completata con processi automatici e richiede un controllo attento e continuo da parte di competenze qualificate Per superare questo osta colo all applicazione su vasta scala di questa procedura realizzando un processo di tempistiche compatibili con le esigenze della pratica clinica indispensabile individuare una procedura quanto pi possibile automatizzata in cui le manipo lazioni manuali siano ridotte al minimo necessario 132 5 1 La selezione delle metodologia meshless applicata Un idea immediata quindi quella di non manipolare minimamente il dataset CT di partenza provvedendo al pi alla sola selezione della regione di interesse per evitare un inutile appesantimento computazionale Vanno ora sottolineate due circostanze importanti i dataset CT offrono dei valori scalari le unit HU in punti disposti secondo una una griglia di punti a spaziatura regolare con ri ferimento all immagine di figura 5 1 la spaziatura risulta costante nel piano xy avendo le immagini tutte la stessa risoluzione Ax Ay e pu risultare variabile lungo l asse longitudinale di scansione z potendo variare la distanz
185. o della ricerca di una metodologia di analisi automatica del comportamento meccanico dei segmenti ossei nell intento di sviluppare una pro cedura di modellazione numerica automatizzata basata direttamente sulle im magini cliniche l eliminazione della fase di creazione della mesh operazione at tualmente ancora semiautomatica consentirebbe un significativo passo in avan ti nella prospettiva delle applicazioni cliniche La verificata equivalenza tra CM e FEM con interpolazione quadratica resta in questo ambito un elemento di carattere accademico non apportando innovazio ne o immediatamente percepibili progressi di carattere metodologico nella filie ra di modellazione basata sui dati diagnostici In assenza di un guadagno in 131 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle termini di precisione prevalgono i vantaggi offerti da un metodo consolidato e per cui sono disponibili svariati strumenti di calcolo avanzati collaudati e affida bili Gli approcci meshless propongono invece un salto qualitativo nel processo di modellizzazione rendendo idealmente possibile una riduzione del tempo di ana lisi e un ulteriore automazione del processo verificate ovviamente che siano ac curatezza generalit e robustezza altri requisiti essenziali di cui si discusso nel capitolo 3 In questo capitolo verr quindi descritta l applicazione di una metodologia me shless del CM allo
186. o inalterata la posizione AQAA 103 4 L approccio meshless con il Metodo delle Celle mutua tra i vari punti della costellazione quello che cambiato nei due casi la posizione relativa tra sistema di riferimento e costellazione di punti Questa osservazione suggerisce un metodo per affrontare tutti quei casi in cui la singolarit della matrice di interpolazione sia dovuta alla particolare disposizio ne tra sistema locale di riferimento e costellazione di punti siccome il flusso di calore attraverso una superficie invariante rispetto al riferimento si pu pensa re di ruotare il riferimento locale attorno al polo in modo tale da rendere massimo il numero di condizionamento della matrice dei coefficienti dell interpolazione Va infine sottolineato come una volta individuata la rotazione del sistema locale si debba provvedere ad esprimere in questo riferimento gli elementi della matrice dei coefficienti dell interpolazione M poich il bilancio tratta grandezze scalari non vi sono come invece sar illustrato a proposito del problema elastico ulte riori operazioni per riportare le grandezze coinvolte nel bilancio nel sistema di partenza 4 3 1 3 Alcuni test numerici Al fine di verificarne l accuratezza l approccio descritto stato applicato alla so luzione della equazione di Laplace V w 0 su un domino piano quadrato di lato unitario centrato nell origine di un sistema di riferimento cartesiano con assi para
187. o pi rapido necessita di una complessa fase di pre processing che richiede un impegno qualificato continuativo fatto che limita intrinsecamente il numero di analisi approntabili da uno steso operatore in un dato tempo Oltre alla validazione del modello MCM l elaborato ha maturato alcuni risultati di carattere teorico nello sviluppo del Metodo delle Celle l estensione ai problemi scalari e vettoriali tridimensionali con interpolazione quadratica ha confermato come il CM possa ottenere ordini di precisione e convergenza paragonabili con 176 quelli ottenibili da sofisticati codici commerciali implementanti il FEM sono poi stati sviluppati e verificati alcuni approcci meshless per la soluzione di problemi bidimensionali scalari e vettoriali che conservando i vantaggi dell assenza del la creazione di una mesh superano alcuni limiti che caratterizzano molte delle formulazioni analoghe in ambito differenziale 177 178 Bibliografia 1 Mark S Shephard Update to Approaches to the automatic generation and control of finite element meshes Appl Mech Rev ASME 49 10 S5 S14 October 1996 2 P L George Automatic Mesh Generation Application to Finite Element Methods John Wiley amp Sons Inc New York NY USA 1991 3 D Testi M Viceconti A Cappello and S Gnudi Prediction of hip fracture can be significantly improved by a single biomedical indicator Ann Biomed Eng 30 6 801 7 2002 4 M C Hobatho R Darm
188. ogni modello pari a 3 mm nella diafisi e di 2 mm nell epifisi il numero di elementi degli otto modelli FEM varia tra 69272 e 80508 il corrispondente numero di nodi tra 104020 e 124592 La distribuzione disomogenea delle proriet meccaniche stata mappata su gli elementi dei modelli FEM tramite il software BoneMat_V3 in grado di attri buire un valore medio del modulo di elasticit E ad ogni elemento della mesh l algoritmo converte dapprima le unit HU in valori di E e poi calcola un integra zione numerica entro il volume spazialmente registrato con l acquisizione CT di ciascun elemento del modello FEM 6 1 3 Acquisizione CT registrazione spaziale e propriet mecca niche I femori sono stati sottoposti a tomografia computerizzata scanner HiSpeed General Electric Co USA parametri di acquisizione 120 kVp 160 mA tipici di un normale esame clinico dopo essere stati scongelati ed immersi in acqua per prevenire artefatti di indurimento del fascio di radiazioni stato applicato un protocollo di acquisizione elicoidale con passo di avanzamento pari a 1 3 mm dalla testa al piccolo trocantere e di 1 5 mm per la regione diafisaria e spessore del fascio rispettivamente pari a 1 e a 5 mm L elaborazione ha restituito un dataset ricostruito con slices di passo 1 3 mm nella zona tra la testa e il piccolo trocantere e di 5 mm nella regione diafisaria e dimensione del pixel pari a 0 59 mm L equazione di calibrazione densitometric
189. olazione locale con funzioni di base radiale essere S 2td 0 da cui d Qi r gu d y dy d 7 3 4 76 calcolando i termini gt Qi2 r Gy x d dx d d Qis kt ge d y dy 4 77 d d Qis ut g t d dx d conformemente all orientazione delle singole superfici Sommando i contributi di flusso calcolati per tutte e quattro le superfici si ottiene Qi aar Qia Quo Lig Qia is Pe 6y 6a 4d 3 Wi 4 78 AiP Wi Il termine appena calcolato esprime il contributo della funzione radiale di centro i al flusso di calore Qaq attraverso la superficie di contorno della regione tribu taria Qp del nodo P assunto come polo Considerando i contributi delle funzioni radiali di centro tutti i nodi della costellazione C si ha pertanto Wi Qanp lap Q2P an p E apWc 4 79 WN e dunque ricordando la 4 69 Qoap ap Te 4 80 che esprime il flusso netto attraverso la superficie di contorno della regione tri butaria del nodo polo in funzione delle temperature dei nodi che compongono la costellazione Il vettore kp ap 4 81 119 4 L approccio meshless con il Metodo delle Celle hy Figura 4 21 Regione tributaria quadrata attorno al generico nodo interno P contiene i valori da collocare nella matrice fondamentale K del sistema che espri me la soluzione del problema KT Q 4 82 alla riga corrispondente all identificativo del nodo po
190. olazione polinomiale con le uniche differenze di carattere puramente implementativo legate alla diversa costituzione della funzione interpolante 4 4 1 Interpolazione RBF di campi scalari In questo paragrafo verr presentata l applicazione della interpolazione con le Funzioni di Base Radiale RBF nell approccio meshless del Metodo delle Celle applicato ai problemi di conduzione termica in regime stazionario La procedura come precisato in precedenza assolutamente simile a quella seguita a proposi to dell interpolazione polinomiale Assunto un dominio di geometria e propriet conduttivit termica k note si provveda a posizionarvi all interno un conveniente numero di punti di seguito denominati nodi La posizione di questi punti verr riferita ad un sistema di riferimento cartesiano esterno XOY Preso in esame un nodo interno generico chiamato convenzionalmente nodo polo si selezionino Nc nodi prossimi disposti attorno al nodo polo questi nodi selezionati prenderanno il nome di nodi satellite e assieme al nodo polo costi tuiranno la costellazione locale C di punti con cui verr approssimato localmente il campo delle temperature attorno al nodo polo 4 4 1 1 Interpolazione Nella regione circostante il nodo polo il campo delle temperature viene appros simato interpolando i valori nodali di temperatura della costellazione C con la funzione Ne T P e X wie o P Pill 4 66 i 1 Per co
191. om pletata da una digressione sulle applicazioni nello specifico ambito della biomec canica I caratteri generali di un approccio meshless reso accessibile in modo diretto dalla struttura concettuale del CM vengono poi declinati nello specifico di tre approcci proposti il primo basato sulla creazione di un complesso locale di celle primali gli altri due basati su un interpolazione della funzione incognita definita localmente su una costellazione di punti a connettivit non definita Pre sentata la formulazione teorica e alcune considerazioni sugli aspetti implementa tivi gli approcci proposti vengono applicati a problemi di test di riferimento sia relativamente alla soluzione nuemrica dell equazione di Laplace che alla soluzio ne di problemi di elasticit lineare nel caso assunto come riferimento della Trave di Timoshenko Nel capitolo 5 viene presentato lo studio di applicabilit di uno degli approcci proposti il metodo selezionato quello basato sulla definizione di un comples so locale di celle primali e viene denominato Meshless Cell Method MCM La 4 verifica di applicabilita avviene replicando con il modello MCM una prova spe rimentale di carico su un campione di femore in vitro problema usato come ri ferimento per la sua valenza rispetto alla prospettiva di applicazione in ambito clinico La valutazione viene condotta comparando le misure sperimentali di de formazione con i valori predetti dal modello numerico ne
192. omogeneo e isotropo lI v 0 E v 1 0 D Eo Eoy 4 9 7cfr 61 90 4 2 Meshless a celle locali con il CM con E modulo di elasticit e v modulo di Poisson del materiale B la matrice che esprime il vettore gradiente degli spostamenti nella cella c 1 m 0 Aa 0 Aa 0 Co OF Ae 0 Ahy 0 Ay 0 Ajy 4 10 Aha Ahy Aix Aig Age Ajy con t spessore del sistema piano A area di cella e in cui compaiono le componenti dei vettori area di cella di cui alla notazione di figura 4 9 Un Un Vi le componenti dello spostamento nodale ai nodi di cella h ej e si scrive il bilancio delle forze al nodo polo imponendo la condizione di equilibrio tra le forze di volume Fp agenti sulla sua regione tributaria e le forze agenti sulle sue superfici di contorno avendo calcolato cella per cella del complesso locale tutti i termini delle forze di superficie T attraverso le porzioni di superficie A di contorno della regione tributaria del nodo polo P siha So Te Fe 4 11 cET P dove il termine a primo membro espresso in funzione dei valori nodali delle componenti di spostamento di tutti i nodi satellite S della costellazio ne Si ottiene una espressione nella forma kea 0 thee O Weg 0 a bier 0 08 F O hee O Kau dg De dop kesal 4 12 I coefficienti della matrice a primo membro andranno a posizionarsi nel la matrice fondamentale del sistema lineare risolvente alle righe associate al nodo polo e ne
193. one del complesso duale ogni cella primale viene ad essere divisa in 10 parti una per ogni nodo Ogni parte di fatto l intersezione fra la cella pri male in esame e le celle duali associate ai suoi nodi Ogni parte rappresenta cos la regione tributaria di ogni nodo nella cella considerata Per scrivere ad esempio il bilancio delle forze agenti sul nodo h si considerano dunque le forze agenti sulla sua regione tributaria Le forze agenti sono e le forze di volume tipicamente la forza peso B nella porzione raccolta dal volume di cella associato al nodo ossia la parte di cella duale del nodo h che cade nella cella c in esame e le forze di superficie Te che agiscono sulle superfici di confine fra la re gione tributaria del nodo A e le regioni tributarie dei tre nodi adiacenti nel caso del nodo h per la convenziona adottata sulla denominazione locale dei nodi i nodi adiacenti sono i nodil men e le forze concentrate F eventualmente applicate direttamente al nodo h o che cadono all interno della sua regione tributaria si noti che per quanto concentrata possa essere l applicazione di una forza essa agisce comunque su una regione che per quanto piccola resta pur sempre di estensione finita il concetto di forza concentrata in un punto che dal punto di vista fisico una pura astrazione logica resta tale anche nell approccio numerico del Metodo delle Celle che attribuisce a buon diritto le forze concentrate alla cella duale
194. opriet materiali avvie ne direttamente dalle informazioni del dataset e l imposizione delle condizioni di vincolo e di carico risulta immediata L applicabilita di questo approccio denominato Meshless Cell Method MCM stata verificata replicando una prova sperimentale di carico su un segmento di femore in vitro Il modello stato costruito su un dataset CT sottocampionato per ragioni di risorse computazionali disponibili del campione osseo registrato rispetto al riferimento sperimentale i valori principali di deformazione registrati in 13 siti anatomici significativi in 5 distinte configurazioni di carico in regime elastico sono state lo strumento di valutazione Confrontando i valori ottenuti con il modello MCM nei punti corrispondenti alla posizione degli strumenti si verificata una elevata correlazione R 0 85 delle predizioni il coefficiente angolare della retta di regressione valori numerici valori sperimentali risul tato pari a 1 21 con intercetta pari a 36 pe lo scarto quadratico medio rispetto alle misure sperimentali risultato pari al 16 con un valore massimo pari a 1 35 5 Verificatane l applicabilit la modellazione MCM stata validata rispetto alle misure sperimentali di deformazione di prove meccaniche in campo elastico con dotte su 4 paia di femori in vitro Misure di deformazione in 16 distinti siti anato mici e in 6 differenti configurazioni di sollecitazione hanno costituito l elemento di valu
195. ore medio delle distanza dei nodi satellite dal no do polo qui conveniente adottare un riferimento polare locale centrato nel no do polo avente asse di riferimento l asse delle ascisse del riferimento cartesiano locale come illustrato dalla figura 4 13 Essendo 0 deo 4 22 y psin il flusso netto attraverso 00 p assume la forma 2T Qoap uf 2dcos R b gy Rcos0 2e sind R c Rsin0 d 4 23 0 2rtk e d R essendo t lo spessore del sistema piano Nodi di contorno Se il nodo polo P risulta posizionato sul contorno del domi nio la procedura descritta resta inalterata Con riferimento alla figura 4 14 e il flusso di calore uscente dalla superficie di confine della regione tributaria assumer in termini generici la forma 02 Qoap te 2 dcos R b gz Recos 2e sin R c Rsin6 d0 01 4 24 e l equazione di bilancio contemplera anche la quota parte dell eventuale flus so imposto al contorno intercettata dalla regione tributaria Qaar Qap Qgnap 4 25 Se il flusso unitario entrante imposto al contorno si ha Osio 2R 4 26 4 3 1 1 Aspetti implementativi Nell ottica di una implementazione della procedura in un codice numerico dedi cato conviene accennare all impostazione vettoriale della procedura descritta 99 4 L approccio meshless con il Metodo delle Celle Figura 4 14 Regione tributaria per un polo situato sul bordo del dominio
196. ore uscente dalla superficie OQ p di confine della regione tributaria 2p del nodo polo P elaborando la funzione che esprime il campo locale delle temperature attorno al polo poich gli operatori matematici che intervengono in questo processo derivazione inte grazione e operazioni algebriche sono lineari ne consegue che l intero processo che calcola il flusso netto attraverso OQp dalla funzione T x y anch esso li neare Se dunque la funzione T x y viene espressa come somma di funzioni componenti 7 il flusso totale Q pu venire calcolato sommando i valori Q otte nuti applicano il processo ad ogni singola funzione addenda 7 Alla luce di questa considerazione tratteremo di seguito una sola delle funzioni che compongono la funzione approssimante corrispondente alla funzione T x y wic Ery Wiec 1 x xi y yi 4 70 centrata nel nodo generico i della costellazione C Le componenti del vettore gradiente della temperatura gc in un punto P x y sono esprimibili in forma analitica derivando la funzione di interpolazione locale 21i pedice c si rende necessario riferendosi all interpolazione locale della costellazione C 117 4 L approccio meshless con il Metodo delle Celle Ti x y ge SU L i una Ea ee da 4 71 Gy on ey 4 1 Sy Yi x y yi wie y Il flusso unitario di calore q correlato con il vettore gradiente delle temperature a meno della conducibilit
197. ose complesse e delicate oltre che 7 1 La creazione di modelli numerici subject specific di segmenti ossei in genere costose prove sperimentali e poter studiare gli effetti di differenti con dizioni di carico su un medesimo sistema La costruzione di un modello subject specific possibile partendo dalle infor mazioni disponibili dalle immagini di una tomografia computerizzata CT tec nica di diagnostica per immagini che allo stato attuale costituisce la migliore e pi attendibile fonte di informazioni circa la morfologia e le propriet meccaniche dei segmenti scheletrici di un soggetto in vita I dati ottenuti con un acquisizione CT consentono di ottenere una rappresen tazione geometrica e in prima approssimazione delle propriet meccaniche di un osso cui pensabile l applicazione di un qualche metodo numerico per l a nalisi delle sollecitazioni fissati opportunamente carichi e condizioni di vincolo Questa operazione comunque tutt altro che semplice innanzitutto per le diffi colt che si incontrano nella costruzione di un modello geometrico in secondo luogo per la complicata caratterizzazione del comportamento meccanico del ma teriale e quindi in ultima analisi per le difficolt nell individuare un modello di comportamento attendibile Molte sono le strade possibili e molti sono i metodi ad oggi sperimentati l inda gine sulla validit di un metodo si basa evidentemente sul confronto con i dati
198. otropo gli elementi della matrice di conduttivit termina K divengono tra loro identici e l equazione 3 27 assume la forma de kgx qy kgy qz kg 3 29 o in forma vettoriale de k Sc 3 30 nota come equazione di Fourier 3 1 1 5 Flusso di calore Il vettore gradiente g e il flusso di calore q sono funzioni affini delle coordina te di conseguenza il flusso di calore Q attraverso una superficie A esprimibi le come prodotto scalare del vettore area che definisce la superficie per il vet tore flusso unitario q valutato nel baricentro G della superficie stessa Essendo A A A A il vettore area e G il suo baricentro si ha Q A Aq A Ay Az 4 G gt 3 31 3 1 1 6 Bilancio L equazione fondamentale del problema viene ottenuta scrivendo il bilancio ter mico riferito alle regioni tributarie associate ad ogni nodo del sistema Queste regioni tributarie sono offerte dalle celle duali Siccome le regioni tributarie sono 53 3 Formulazione quadratica con il Metodo delle Celle arbitrarie per forma e per estensione anche il criterio di definizione delle celle duali anch esso arbitrario Questo significa in termini pratici che tutti i crite ri di divisione delle celle tetraedriche in porzioni associate ai suoi nodi sono di principio ammissibili ancorch almeno dal punto di vista computazionale non equivalenti Nel caso bidimensionale con i simplessi ad interpolazione quadra
199. rare la validit del criterio adottato e la manipolazione manuale del dataset ha evidenziato come il valore di HU pari a 68 non escluda porzioni di tessuto osseo anche nelle regioni a minore densit media quali quelle prossimali in cui l osso ha struttura trabecolare e la variazione del valore di soglia nell intorno del valore adottato ha eviden ziato una sostanziale stazionariet dei valori di deformazione calcolati in sale calcio fosfato che rappresenta la frazione minerale tessuto osseo dove presente in forma cristallina amorfa o mista 159 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle corrispondenza degli estensimetri il grafico di figura 6 7 evidenzia come al variare del valore HU i valori di deformazione calcolati relativamente ad uno dei femori in una stessa condizione di sollecitazione si mantengano sostanzialmente stabili nell intorno di HU 68 6 2 2 L individuazione dei valori puntuali di deformazione Sia con che senza soglia HU nell attribuzione del modulo di elasticit i valori predetti in alcuni estensimetri sono risultati sistematicamente molto superiori sia alle misure sperimentali che alle corrispondenti predizioni dei modelli FEM A fronte di una base di valori confrontabili queste misure risultavano sistematica mente anomale sia nei valori in s che nel comportamento globale del modello MCM Verificando il posizionamento di questi punti sui datase
200. re in luogo del riferimento cartesiano globale un riferimento locale per ciascuna cella Un riferimento locale offerto dalle coordinate affini locali definito dai tre lati uscenti da uno stesso vertice 48 3 1 Problemi scalari tridimensionali Figura 3 1 Cella tetraedrica a 10 nodi disposizione e denominazione convenzio nale dei vertici e sistema di riferimento locale di coordinate affini 3 1 1 1 Riferimento locale di coordinate affini Il legame tra le coordinate locali affini e le coordinate cartesiane globali lineare x Li j i k xi ci x E YP Suip y vu M np 3 2 zZ Zi i zi Zk _ zi z Zi che con una sintetica notazione vettoriale di evidente significato pu venire espressa come r r J s 3 3 La relazione inversa dunque r a B y fz ng 4s e O lt y 3 4 t UVT z 49 3 Formulazione quadratica con il Metodo delle Celle essendo a 1 Allu v a zi yw ys Ze zi B 1 A ai Ti Ze 21 k x z1 2 y 1 A p 2 yi yi t Fi Ye y 6 1 Al y yi z de yj yi a 2a e 1 A x xi a a1 vi z 9 1 A a x yj ys c 2 w va pe Ally yi 2r 2 Yk Yi 5 2 v IA eem m 2 21 y 2 a m 1 A z xi Yk yi ar y vl e Ta CU Lk Y
201. re e al diametro della testa L assenza di un estensimetro nell aspetto laterale del collo dovuta all irregolarita della superficie in questa regione 151 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle 6 1 1 1 Condizioni di carico Su ciascun femore sono state applicate sei differenti condizioni di carico corri spondenti a flessione in piani differenti carico assiale e torsione 4 configurazioni di carico corrispondono alle direzioni estreme raggiunte dalla forza di reazione all anca nello svolgimento di normali azioni motorie 79 ulteriori configurazioni corrispondono al carico verticale direzione definita dal piano sagittale da quello frontale del femore e ad una condizione inclinazione del femore di 8 nel pia no frontale assimilabile alla fase di appoggio monopodalico della camminata e individuata essere con una preliminare indagine FEM particolarmente onerosa per la regione del collo femorale in termini di stress prodotti Al fine di prevenire danni meccanici ai campioni stato applicato un carico mas Figura 6 2 Visualiz Figura 6 3 Mesh di un modello FEM zazione delle varie nella regione prossimale con evidenziati configurazioni di ca gli elementi afferenti ai nodi corrispon rico denti ad estensimetri simo corrispondente al 75 del peso corporeo del donatore Ciascun campione stato montato sulla cella di carico di portata 5k N della macchina di prova Mod 8502
202. riscono un approccio locale alla costruzione della equazione di bilan cio Dal punto di vista del singolo nodo primale preso in esame sulla cui circo stante regione tributaria cella duale si sta scrivendo il bilancio la conoscenza dei due complessi di celle nella loro globalit evidentemente superflua dal punto di vista del nodo esaminato infatti sufficiente conoscere la connettivit locale tra questo nodo ed i nodi primali circostanti connettivit che permette di definire l insieme delle celle primali che afferiscono al nodo stesso Definite che siano localmente queste celle infatti possibile definire sempre localmente gli elementi del complesso duale ad esse associati questo consente in particolare di definire la cella duale che circonda il nodo considerato ossia la sua regione tributaria rispetto a cui esprimere il bilancio di nodo Questa procedura in cui la divisione delle stesse celle primali in porzioni duali avviene pi volte per cia scuno dei vertici di cella ridondante e computazionalmente non ottimizzata se applicata alle celle di un complesso definito nella sua globalit a priori Il pro cesso solitamente seguito quello di esaminare una cella alla volta calcolando porzioni del bilancio di nodo dei suoi vertici conducendo questa operazione per tutte le celle primali si coinvolgono tutte le porzioni di cella duale di competenza di ciascun nodo e il loro corretto assemblaggio porta immediatamente
203. ro satelliti Figura 4 22 Interplazione meshless con RBF errore quadratico medio al variare del numero dei satelliti con le varie funzioni di base radiale e passo di griglia pari a 1 8 quelle proprie delle altre tre funzioni e stabili rispetto all aggiunta di un ulteriore nodo nella costellazione il numero dei satelliti deve essere di almeno 8 I dati di tabella 4 8 riportati nel grafico di figura 4 23 riportano invece il valore di 6 al variare del passo p della griglia di punti nelle loro disposizione regolare di origine lato 6 a bs oc 0 p 0 5 0 009267 0 006156 0 004561 0 004013 0 25 0 002386 0 001551 0 000815 0 000603 0 125 0 000583 0 000400 0 000175 0 000118 0 0625 0 000143 0 000096 0 000034 0 000022 0 03125 0 000035 0 000024 0 000010 0 000001 Tabella 4 8 Meshless con interpolazione RBF scarto quadratico medio 6 al variare del passo griglia per le varie funzioni di base radiale con costel lazione composta da 8 satelliti 122 4 4 Interpolazione locale con funzioni di base radiale 24 e n e n 34 A e v n A v 44 n i D 2 e A n v 54 A oA E oB A oC vod 6 v T T T T T T 1 1 6 1 4 1 2 1 0 8 0 6 0 4 0 2 log lato Figura 4 23 Meshless con interpolazione RBF scarto quadratico medio al variare del lato in scala logaritmica per le quattro funzioni di base radiale 4 4 2 Interpolazione RBF di campi vettoriali In questo paragrafo verr present
204. rodotte dalla misura dalla trasmissione dei raggi di X X ray Computed Tomogra phy Tuttavia il termine potrebbe anche includere le immagini ottenute con la SPECT Single Photon Emission Computed Tomography o con la PET Positron Emission Tomography e perfino alcune tecniche ad ultrasuoni basate su misure di proiezione Le immagini CT sono tipicamente degli array bidimensionali di 512 512 pixel Ogni pixel picture element elemento bidimensionale di base di un immagine digi tale riporta un valore scalare che individua uno dei 4096 possibili gradi di grigio definibili con una codifica a 12 bit Siccome le immagini slices sono tra loro sepa rate ad ogni immagine associato uno spessore thickness dell oggetto scandito e dunque ad ogni pixel associato un voxel volume element Il vantaggio principale di una CT rispetto alle proiezioni semplici la capacit di separare gli oggetti disposti nel senso della proiezione evitando la confusio ne che si crea quando le ombre di pi oggetti si sovrappongono Questo fatto unitamente all alta precisione delle misure e alla produzione di immagini digita li d alla CT la capacit di cogliere oggetti con contrasto estremamente piccolo Per esempio i sistemi CT che usano dei raggi X sono in grado di distinguere fa cilmente oggetti per i quali la differenza relativa nel coefficiente di attenuazione una frazione di per cento Questo permette rappresentazioni altrimenti non ottenibili con le
205. rtendo da dati incompleti affetti da rumore L accuratezza di in un metodo pu essere valutata in vitro ma allo stato attuale difficilmente accertabile in vivo Di fatto esiste la possibilit di misure in vivo di micro motion per mezzo della Rontgen Stereophotogrammetric Analysis questa tecnica stata usata principalmente per la valutazione di impianti micro motion ma la sua accuratezza ritenuta discutibile ed si ancora nella fase di standar dizzazione del metodo Il metodo di calcolo attualmente pi affermato senza dubbio il Metodo de gli Elementi Finiti FEM gi ampiamente utilizzato in ambito biomedico per lo studio di problemi di campo termico magnetico elettrico e per l analisi della diffusione delle radiazioni nei tessuti corporei Un altra applicazione quella ri guardante problemi strutturali per stabilire le deformazioni e gli stati di tensione indotti da azioni esterne sugli organi del corpo umano ivi compreso pertanto l ambito ortopedico Il FEM si affermato rispetto ai metodi concorrenti Boundary Elements Me thod Finite Volume Method per la sua attitudine a studiare sistemi di ele vata complessit geometrica e meccanicamente non omogenei eventualmente caratterizzate da un comportamento non lineare 1 3 La definizione di un modello In letteratura vengono riportate molte metodologie per la definizione di modelli ad elementi finiti subject specific caratterizzate da diversi liv
206. rticolarmente nella zona del col lo a causa dell inclinazione relativa in tale zona della superficie e del fascio di 141 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle Figura 5 9 Modello virtuale tridimensionale NURBS del femore radiazione incidente La regione condilare non stata inclusa nel modello FEM perch non si ritenuto necessario modellare la zona inclusa in PMMA stiman do sufficiente il vincolo fisso dei nodi giacenti al livello prossimale del blocco La mesh ottenuta risultata di 76026 elementi e 118970 nodi La mappatura delle propriet meccaniche del tessuto osseo sugli elementi del la mesh avvenuta con il software BoneMat 81 Il software attribuisce ad ogni elemento della mesh un valore di modulo di elasticit sulla base dei valori di uni t HU lette nel dataset CT nella regione di spazio occupata dall elemento Oltre al dataset CT il software richiede come input l equazione di calibrazione densi tometrica legame HU densit dell acquisizione CT e la relazione costitutiva legame densit modulo elastico entrambe descritte nel seguente paragrafo 5 2 5 Il modello FEM risulta quindi elastico lineare isotropo nel comportamen to materiale di ciascun elemento ma disomogeneo nella distribuzione dei moduli di elasticit 142 5 2 Il test esplorativo 5 2 3 1 Condizioni al contorno Per riprodurre le condizioni di vincolo d
207. screto non produce singolarit formali nell applicazione di sorgenti concentra te n sussistono condizionamenti nella forma del dominio analizzabile o vincoli nella discretizzazione Esso si candida dunque ad essere un metodo potenzial mente utile per lo studio di sistemi biomeccanici tipicamente caratterizzati da geometrie complesse e irregolari e da una grande variabilit distribuzione delle propriet materiali L esplorazione di alternative al FEM stata focalizzata sul Metodo delle Celle per due ragioni in primo luogo il CM rifiutando a priori l ipotesi di continuo e assu mendo filosoficamente la prospettiva dell approccio discreto appare nei principi concorde con la struttura della fonte di informazione i dataset di immagini dia gnostiche con cui i modelli possono venire definiti in secondo luogo il CM ha manifestato di poter rendere accessibili alcuni strumenti quali la formulazione quadratica con ordine di convergenza nativamente pari a 4 e un approccio me shless di semplice formulazione tali da costituire a priori motivo di interesse e ragione di studio in un ottica non gi solamente accademica quanto piuttosto pra tica e applicativa Dal punto di vista concettuale la tesi strutturata in tre parti la prima par te inquadra il problema di riferimento e presenta il metodo numerico oggetto 3 di studio la seconda parte focalizzata sullo sviluppo formale di alcune imple mentazioni innovative la terza parte ha
208. si le equazioni ottenute con l approccio differenziale descrivono il comportamento dei sistemi analizzati in modo esatto rispetto alle ipotesi poste questa soluzione esatta per quantitativamente accessibile solo in pochi e particolari casi per geometrie in genere regolari ed in presenza di parti colari condizioni al contorno Per poter analizzare sistemi reali le equazioni differenziali devono essere di scretizzate fatte risalire cio dall infinitesimo a dimensioni che pur convenien temente piccole sono comunque finite Si ottengono in questo modo dei sistemi di equazioni algebriche oggi facilmente risolvibili con l utilizzo del calcolatore elettronico Per questo scopo sono state elaborate varie metodologie alcune tra le 25 2 Il Metodo delle Celle pi note delle quali sono il Metodo delle Differenze Finite FDT il Metodo degli Elementi Finiti FEM il Metodo degli Elementi al Contorno BEM e il Metodo dei Volumi Finiti FVM L implementazione di queste procedure in genere complessa ed ha il difetto di perdere il legame con il problema fisico in esame Il Metodo delle Celle segue invece una strategia differente anzich affrontare le equazioni differenziali fondate sull astrazione teorica del punto per descrive re una realt fisica discreta il Metodo si basa sulla formulazione discreta diretta delle leggi fisiche Fissando l attenzione direttamente su porzioni di spazio finite le leggi fi
209. siche del fenomeno in esame sono scritte cio direttamente in termini discreti partendo dalle leggi sperimentali Questa strada permette di costruire direttamente in forma algebrica le equazioni del campo fisico studiato equazioni di bilancio equazioni circuitali equazioni materiali Il metodo conserva l ade renza alla realt fisica del problema studiato e al contempo segue la metodologia dell approccio sperimentale sia nella scala di osservazione che nell attribuire le grandezze fisiche agli elementi geometrici dello spazio Questi argomenti meritano una specifica digressione poich costituiscono la pre messa concettuale della formulazione del Metodo delle Celle 2 1 Le variabili fisiche Lo studio delle grandezze fisiche la prima e fondamentale fase che accompagna la formulazione di un modello descrittivo di qualunque fenomeno fisico Vo lendo inquadrare la formulazione finita di una legge fisica importante partire proprio da una riflessione attorno alle grandezze fisiche alla loro natura e alla loro classificazione 2 1 1 Costantie variabili La grandezze fisiche possono dapprima essere divise in due grandi famiglie un primo gruppo di grandezze costituito da quei valori che assumono il ruolo di parametri del problema descrivendone la natura e qualificando attributi il cui valore assunto invariabile nel fenomeno studiato questa la categoria del le costanti fisiche cui appartengono tanto le costanti universali entit
210. solidata applicazione per il quale disponibile una variet di strumenti di gestione e di calcolo evoluti e affidabili Il secondo ambito di indagine ha riguardato lo sviluppo di implementazioni me shless del CM approcci nei quali la soluzione viene calcolata in punti tra i quali non richiesta la definizione di una connettivit a livello globale Sono stati svi luppati nella loro impostazione formale e implementativa tre differenti approcci questi approcci hanno in comune il fatto di scrivere l equazione di bilancio con riferimento a regioni tributarie definite localmente nodo per nodo La prima del le ipotesi proposte usa un complesso locale di celle primali definito attorno ad ogni nodo esaminato le celle locali sono usate sia per definire la regione tribu taria che per il calcolo cella per cella dei termini di flusso coinvolti nel bilancio di nodo Negli altri due approcci la regione tributaria assunta come una figura geometriche di forma elementare il campo incognito approssimato nell intorno di ogni nodo tramite una funzione analitica che interpola i valori incogniti del nodo in esame e di un gruppo di nodi circostanti due ipotesi sono state esplorate circa la funzione interpolante i polinomi ordinati di grado parziale uniforme e le funzioni di base radiale Le formulazioni di questi approcci sono state verificate con dei test numerici di confronto con le soluzioni teoriche sia relativamente ai campi scalari soluzione de
211. sperimentali ottenuti in laboratorio su campioni espiantati 1 2 Criteri di applicabilit In linea di principio l analisi numerica delle sollecitazioni meccaniche in una struttura ossea possibile con un qualsiasi metodo di calcolo Affinch un meto do di analisi numerica possa essere utilizzabile per studi clinici esso deve per soddisfare dei requisiti di base e dovendo risultare e Automatico Automation per ridurre il tempo di analisi il numero del le prove da effettuare in genere elevato e le applicazioni diagnostiche richiedono risposte in breve tempo e Generale Generality per poter rappresentare correttamente ogni segmento osseo le ossa dello scheletro umano presentano infatti una grande variet di forma dimensioni e propriet fisiche e Accurato Accuracy perch i risultati della analisi possano poter essere uti lizzati in fase di diagnosi o di definizione di terapie e Robusto Robustness per riuscire a produrre un modello accurato anche da dati incompleti o affetti da rumore 8 1 3 La definizione di un modello Per meglio precisare questi concetti Automaticita attitudine a completare il processo di analisi con il minor numero possibile di interventi Generalit attitudine a completare il processo di analisi per qualunque seg mento di osso a prescindere dall eventuale presenza di anomalie anatomi che fisiologiche o patologiche Robustezza attitudine a condurre analisi anche pa
212. spondente ad un suo ottavo ottenuto sezionando la sfera con i piani coordinati come illustrato in figura 3 8 Per replicare le condizioni cinematiche del problema originario sui punti del le superfici di taglio si sono eliminati gli spostamenti ortogonali alla superficie stessa consentendo solamente spostamenti nel piano della sezione per i punti appartenenti agli spigoli radiali gli spostamenti vincolati saranno pertanto risul tati dall unione dei vincoli imposti sulle superfici concorrenti Questa geometria stata discretizzata con elementi tetraedrici a 10 nodi otte nendo 3 differenti mesh a differente dimensione caratteristica l ottenuta come rapporto tra la lunghezza di uno spigolo radiale e il numero di divisioni su di esso impostate La figura 3 9 illustra una delle discretizzazioni utilizzate nel l analisi Lo stesso complesso di elementi stato utilizzato per risolvere lo stesso problema con il programma Ansys codice commerciale general purpose per la soluzione numerica di problemi ingegneristici basato sul FEM Per ciascuna delle discretizzazioni si valutato in ogni nodo non vincolato il valore dello spostamento radiale u confrontandolo con il valore teorico i la regolarit della geometria e l attitudine della tassellizzazione di Delaunay alla base del discretizzatore a campire una geometria di questo genere hanno fatto s che le mesh fossero rego lari e di buona qualita in queste condizioni la
213. sso usa le funzioni di base radiale per l interpo lazione locale delle grandezze incognite e viene qui presentato e descritto re lativamente ad entrambi i problemi sin qui trattati campi scalari e problema dell elasticita piana Una prima e immediata idea per interpolare localmente nell intorno del nodo polo i valori assunti dalla variabile di configurazione incognita nei nodi circo stanti l adozione di una funzione polinomiale Ragioni di simmetria il proble ma prescinde dal riferimento adottato impongono la scelta di un polinomio in cui le coordinate coinvolte figurino con uguale grado Specificamente la scelta avviene attingendo al cosiddetto triangolo di Pascal per la selezione delle fun zioni interpolanti Nel seguito verranno descritti nel dettaglio le formulazioni per i due problemi scalare e vettoriale di riferimento 4 3 1 Campo scalare stazionario Considerando come riferimento il problema della conduzione del calore in regi me stazionario all interno di un dominio bidimensionale sia dato un sistema pia no nel quale siano stati posizionati un certo numero di punti di seguito chiamati nodi la cui posizione sia riferita ad un sistema di assi cartesiani di riferimento Lo scopo quello di determinare la temperatura in ciascuno dei punti note che siano le condizioni al contorno temperature o flussi di calore imposti Si consideri un generico nodo P interno al dominio e si adotti un sistema di ri feriment
214. t CT si rilevato come essi fossero collocati subito all esterno dell osso con errore di qualche pixel nel l immagine o in regioni di transizione tra tessuto osseo e matrice esterna acqua che per la risoluzione spaziale della tecnica di acquisizione delle immagini ri sultano zone con densit a cavallo tra quella dell osso e quella dell acqua Questo fatto direttamente imputabile alla combinazione di due errori nel posi zionamento dei punti di misura posizione sperimentale degli estensimetri nel riferimento del dataset CT e l incertezza nel rilevamento della posizione degli estensimetri sulla super ficie del femore operazione che avviene con uno strumento digitalizatore di risoluzione 0 2 mm e la registrazione tra il riferimento sperimentale e quello del dataset CT ope razione caratterizzata da una incertezza dell ordine dei 0 5 0 9 mm Il casuale combinarsi di questi errori ha dunque ragionevolmente prodotto ler rata collocazione dei punti di indagine nel modello MCM che avendo la risolu zione nativa del dataset CT risulta estremamente sensibile la risoluzione nativa dell acquisizione CT pari nel piano di ogni singola slice a 0 59 mm Va rilevato come nel modello FEM questo problema non si ponga per due ragioni sostanziali innanzitutto in corrispondenza della posizione stimata dell estensi metro sul femore viene collocato sulla superficie del modello FEM un nodo il quale a prescindere da errori di pos
215. ta ragionevole pensare che un miglioramento delle predizioni del rischio di frattura possa provenire dall uso di modelli numerici tridimensionali modelli subject specific costruiti sulle immagini medicali sono infatti di principio in grado di contemplare in maniera pi completa i fattori determinanti il comportamento meccanico di ogni singolo segmento osseo studiato Requisiti fondamentali di una metodologia di indagine numerica nella specifica prospettiva delle applica zioni cliniche sono oltre all accuratezza la robustezza nei confronti degli artefatti nelle immagini diagnostiche la generalit rispetto al sito di indagine e l automa ticit del processo Il Metodo degli Elementi Finiti FEM a tutt oggi uno dei metodi numerici pi affermati per l indagine e la caratterizzazione del comportamento di un si stema fisico sotto determinate e controllate ipotesi Questo vieppi vero nello specifico ambito della biomeccanica ortopedica in virt della sua attitudine a re plicare accuratamente sia la morfologia che le propriet materiali di un segmento osseo partendo dai dati di una Computed Tomography CT Gli studi svolti nel 1 corso degli negli ultimi vent anni hanno portato allo sviluppo di una metodo logia di modellazione basata sul FEM che partendo da immagini medicali si dimostrata in grado di ottenere precisioni compatibili con le applicazioni cliniche nella replicazione di misure sperimentali di deformazioni su
216. tazione della precisione del metodo nell identificare la realt sperimenta le I modelli MCM di ciascun femore sono stati costruiti sui corrispondenti da taset CT a risoluzione nativa densitometricamente calibrati L analisi aggregata ha confermato la buona capacit dell MCM nel replicare le misure sperimentali il coefficiente di correlazione tra valori predetti e valori sperimentali risultato 175 elevato R 0 94 la retta di regressione con pendenza pari a 0 93 e intercetta pari a 1 88 ue avendosi uno scarto quadratico medio pari del 7 92 e un valore di picco pari al 48 Questi valori sono risultati del tutto paragonabili a quelli ottenuti dalle analisi FEM associate alle misure sperimentali Per comprendere il comportamento del modello MCM rispetto al metodo FEM di riferimento stato approntato un confronto sistematico con uno dei modelli FEM l analisi avvenuta comparando gli spostamenti ed i valori principali di defor mazione ai centroidi degli elementi del modello FEM con i valori calcolati dal modello MCM nei punti corrispondenti L analisi aggregata ha dimostrato una piena corrispondenza tra i due modelli nel dettaglio la concordanza nella defini zione del campo di spostamenti avvalora l ipotesi che la minima dispersione nei valori principali di deformazione sia imputabile alla differente scala di dettaglio con cui i due modelli sono definiti il modello MCM a coglie infatti le prorpiet materiali con la risoluzione n
217. te come il suo segno e quindi anche il suo significato dipenda dall orientazione attribuita alla superficie attraversata Cambiando l orientazione della superficie cambia di conseguenza il segno del flusso ma il suo significato ad esempio di flusso entrante o uscente in un volume di cui la superficie contorno non viene alterato 29 2 Il Metodo delle Celle 2 2 L orientazione degli elementi spaziali Gli elementi geometrici dello spazio utilizzati usualmente per riferire le variabili globali possono tutti venire orientati nello spazio Pi precisamente sono possibi li due tipi di orientazione di un elemento spaziale quella interna e quella esterna ciascuna delle quali dotata di versi 2 2 1 Orientazione interna degli elementi spaziali L orientazione interna di un elemento geometrico un orientazione che viene fat ta sull elemento stesso Da un punto di vista puramente geometrico questo signi fica che per definire un orientazione interna non necessario stabilire l ambiente geometrico in cui l oggetto inserito Si ha cos che un punto viene orientato internamente stabilendo un verso di av vicinamento considerato come positivo ossia stabilendo se il punto sia sor gente oppure pozzo Per orientare internamente una linea se ne stabilisce un verso di percorrenza per un elemento di superficie invece necessario stabilire un verso di percorrenza del suo contorno Un volume viene infine orie
218. te degli spostamenti V u di dv Eyy dy dw Eza 7 Oe 3 45 B E u dU Yay SS ee 7 z So du w _ Ve dz Ox 2 _ B v dw _ a SETT dz y 2 e ee 3 46 che riporta tutte le 6 componenti significative del tensore rappresentato Oltre ad essere questa una notazione comoda essa soprattutto utile ai fini dell im plementazione di un programma di calcolo efficiente In pi nel caso di com portamento elastico lineare del materiale essa permette di esprimere il legame costitutivo con una semplice matrice di dimensione 6x6 anzich con un tensore a quattro indici Chyrs 62 3 2 Problemi di elasticita lineare tridimensionali 3 2 2 2 Gradiente degli spostamenti Il gradiente degli spostamenti viene espresso dalla matrice Or Oy Oz v dv Ov _ PO ed Du 47 Q Or Oy Oz oe dw Ow w Ox Oy Oz Indicando con wu la generica componente del vettore spostamento u e con e la variabile di derivazione si ha evidentemente Qui _ du OE Ou On Ou OC 3 48 de Of deg Onde OC de ER Si ha pertanto du du du 06 dE O 0 amp Om OC x dy dz dv dv Ov On On On arl SS a 2 A Ii 3 49 Q E On OC Oxr Oy dz of dw du w a ac a de On ACI Lox dz Az Va rilevato subito che la seconda matrice altro non che la matrice trasposta del l inversa della matrice J che appare nell equazione 3 3 Analizzando invece la prima matrice appare evidente la scrittura Q ea B EnC u
219. tender appieno attorno al nodo se il nodo giace invece sul contorno del solido in esame la regione tributaria sar la sola parte del cerchio contenuta nel solido Per esprimere le componenti delle forze di superficie e agenti sul contorno della regione tributaria risulta conveniente adottare un riferimento in coordinate polari centrate nel nodo polo e con asse principale parallelo e concorde con l asse delle ascisse del riferimento globale Si ha evidentemente 0 la aed 4 50 y psin 108 4 3 Interpolazione locale con funzioni polinomiali In questo riferimento polare locale le componenti del tensore degli sforzi di vengono A _ E ag sinf R 2a7 cos0 R ag ve 1 v i v E 2a sind R aio cos0 R ag 1 p 4 51 v E ag sind R 2a7 cos0 R az E 2 ai2 sinf R aio cos 0 R ag 1 n i 1 1 Oy 4 52 l v E 2a sin R aio sind R ag cost R 2ag cos R as a4 4 1 v Ti 4 53 La risultante delle forze di superficie agenti sul contorno della regione tribu taria si ottiene integrando opportunamente le tensioni 27 F f R a cos 0 T sin d 4 54 ae bp f R osin 0 T cos d 4 55 0 da cui F 2a11 3 a10 v 8a7 2 a11 a ea 4 56 z 11 10 7 De an F 3 2ag V dg 2ag 8 a12 2 R 4 57 3 ag V y 9 ag ag ag a12 7a eae a Nodo di bordo la regione tributaria del nodo polo deve essere interna al
220. terazioni di posizione 86 4 2 Meshless a celle locali con il CM Figura 4 6 Bilancio sulla regione tributaria di un nodo di bordo quando presente un valore di flusso sul contorno del dominio 0 8 0 6 Figura 4 7 Funzione armonica e cos y nel dominio di lato unitario oggetto del test 87 4 L approccio meshless con il Metodo delle Celle utilizzato il valore dello scarto quadratico medio rispetto alla soluzione teorica Ni 2 S wes wi 4 7 N essendo N il numero di nodi interni quelli su cui non sono stati imposti i valo ri della funzione soluzione del problema w il valore assunto dalla soluzione teorica nel punto i e w il valore ottenuto dalla soluzione numerica nello stesso punto Il problema cos posto stato risolto con la procedura meshless a celle locali so pra descritta implementata con riferimento a costellazioni locali composte da 4 satelliti La stima dell ordine di convergenza ossia del decrescere dell errore al diminuire di una misura significativa della distanza dei punti stata effettuata con riferi mento ad una griglia regolare di punti situazione in cui il passo di griglia d una misura esatta della dimensione associata alla distribuzione di punti La tabella 4 1 riporta i valori ottenuti per il parametro 6 con passi della griglia man mano ridotti lato 0 5 0 002604 0 25 0 000543 0 125 0 000122 0 0625 0 000029 0 03125 0 00
221. teristiche che dipendono dai para metri impostati per la scansione oltre che dalle caratteristiche dello scanner CT Errori di segmentazione La seconda fonte di errore la segmentazione della regione di interesse so no stati proposti vari algoritmi di segmentazione caratterizzati da vari gradi di automazione partendo dall estrazione completamente manuale del contorno e arrivando a complessi algoritmi completamente automatizzati Qualunque algo ritmo di segmentazione si scelga sar in ogni caso presente un incertezza nella definizione del contorno Nella maggior parte dei casi la segmentazione avviene facendo uso di algoritmi basati su processi di sogliatura i cui risultati dipendo no dal valore di soglia impostato 22 6 Questo errore nella definizione della geometria si propaga evidentemente attraverso tutti i passi successivi sino alla creazione della mesh Errori di mappatura In molti studi le informazioni circa la densit ottenute dal data set CT sono utilizzate per derivare una distribuzione non omogenea delle propriet meccani che dei tessuti Questa procedura prevede due passi per prima cosa i dati CT sono calibrati per correlare le Hounsfiel Units HU ai valori di densit apparen te o di ash density del tessuto questa operazione talvolta realizzata utilizzando un fantoccio di calibrazione 23 12 24 5 25 ovvero assumendo dei valori convenzionali attinti dalla letteratura per determinate sezioni selected re
222. ti dello spostamento u Ox Ex Ov Ey sr Oy e 2 dy da 4 47 ag Yy 2a7 z a3 _ 2 a12 y aio ag 2 a11 yY aio y agT 2ag xt as a 2 107 4 L approccio meshless con il Metodo delle Celle 4 3 2 3 Sforzi Per un materiale isotropo a comportamento elastico lineare il legame tra sforzi e deformazioni descritto dalla relazione lineare I v 0 O E if B E v 1 0 i Oy ae 1 Di v Ey 4 48 Try 0 0 5 xy da cui tramite le 4 47 E ag sind R 2a7 cos R az v E 2a sind R aio cos R ag 1 rv i 1 v i _ J vE ag sind R 2a cos R ag E 2a sin 0 R aio cos 0 R ag ae 7 1 v 1 v 1 v E 2a sinf R aio sind R ag cost R 2ag cost R as a4 4 1 v 4 49 4 3 2 4 Bilancio delle forze al nodo Come fatto in precedenza si assume come regione tributaria del nodo polo una circonferenza di centro il polo stesso e di raggio R il valore di R di principio arbitrario e viene qui assunto pari alla met del valore medio della distanza di nodi satellite dal polo 66 Va qui ricordato che in presenza di forze di volume le regioni tributarie devono coprire nel loro insieme l intero dominio studiato in caso contrario laddove restino dei buchi parte dei carichi non venendo raccolta da alcuna regione tributaria e non venendo quindi computata in alcun bilancio resterebbe trascurata Per un nodo interno la regione tributaria si es
223. tica si ve Figura 3 2 Suddivisione del tetraedro in porzioni duali basata sui punti Gauss dei lati rificata una dipendenza tra il criterio di definizione delle celle duali e il grado di convergenza del Metodo In particolare quando la creazione delle celle dua li venisse appoggiata ai punti di Gauss dei lati delle celle primali il metodo ha presentato ordine di convergenza pari a 4 superiore a quanto ordinariamente ottenibile con una corrispondente analisi condotta con il FEM Assunto un criterio per la definizione della suddivisione duale vi sono due possibili modi di procedere e esaminando un nodo alla volta individuare l insieme n delle celle che afferiscono al nodo e definire il poliedro cella duale del nodo che ne costi tuisce la regione tributaria il poliedro duale sar composto da varie porzio ni ciascuna delle quale risulter appartenere ad una cella primale differente 3l equazione di bilancio ha valore a prescindere sia dalla forma che dalle dimensioni della regione a cui riferita 54 3 1 Problemi scalari tridimensionali Figura 3 3 Suddivisione duale per celle triangolari con interpolazione quadrati ca a sinistra la divisione di tipo baricentrico ossia basata sui punti medi dei segmenti congiungenti i nodi a destra la suddivisione si appoggia al baricentro e ai punti di Gauss dei lati di cella Per ogni porzione si calcola il flusso di calore Q A in transito nella cel
224. tiche e dinamiche forze generalizzate forza impulso pressure quantit di moto massa gravitazionale densit di massa momento angolare pressione sforzo entropia flussi termodinamici carica elettrica flusso dielettrico corrente elettrica vettore campo magnetico vettore induz elettrica potenziale magnetico corrente di probabilit ece Classificazione generale delle grandezze fisiche e le variabili di configurazione che definiscono appunto la configurazione del sistema Appartengono alla stessa categoria tutte le variabili ottenute 28 2 1 Le variabili fisiche da variabili di configurazione a mezzo di operazioni di somma differenza moltiplicazione per una costante integrazione derivazione spaziale o tem porale a patto che nelle operazioni non figurino costanti fisiche Nell ambi to della meccanica dei continui sono variabili di configurazione le variabili geometriche e le variabili cinematiche e le variabili di sorgente sono quelle che definiscono le sorgenti di un cam po ossia le entit che generano il campo Per il campo elastico sono variabili di sorgente la massa e le forze Appartengono alla stessa classe tutte le variabili ottenibili da variabili di sorgente a mezzo di operazioni di som ma differenza passaggio al limite derivazione e integrazione divisione per una lunghezza un area un volume e un intervallo di tempo Queste relazioni non devono contenere costanti fisiche e le var
225. to dalla necessit di coprire tutta una regione circostante il nodo polo e dall altro dall esigenza di evitare creazione di celle distorte Zovatto e Nicolini 60 hanno rilevato comunque come l ordine di convergenza ottenuto non dipenda dal numero di nodi satellite utilizzati quando questi varino tra 4 e 8 un numero maggiore di nodi satellite risulta poco conveniente aumen tando la difficolt nella creazione di celle a forma regolare Per ciascuna delle celle del complesso locale si provvede quindi al calcolo del le grandezze coinvolte nel bilancio del vertice corrispondente al nodo polo La procedura cella per cella risulta essere la stessa seguita quando la cella veniva considerata come parte di un complesso globale Con questa procedura la matrice fondamentale del sistema di equazioni algebri che che permettono di risolvere il campo studiato ossia il calcolo dei valori delle sue variabili di configurazione viene costruita riga per riga la riga elaborata quella corrispondente numero associato al nodo polo nella numerazione glo bale dei nodi posizionati all interno del sistema studiato gli elementi non nulli di questa riga saranno quelli posizionati nelle colonne di indice pari al numero identificativo dei vari nodi satellite nella numerazione globale Non appare del tutto superfluo in questo contesto sottolineare come questo ap proccio locale comporti la perdita di garanzie circa la simmetria della matrice fondamentale del sistem
226. tta dal FEM e soluzione teorica al variare del lato di cella della mesh per le due funzioni di test 3 2 1 Interpolazione quadratica del campo degli spostamenti Ciascuna delle tre componenti dello spostamento rispetto ad un riferimento car tesiano globale viene approssimata all interno di ogni cella con il polinomio di secondo grado ui 2 y 2 a bix ciy diz eizy fiyz gizx hz i Liz 3 36 59 3 Formulazione quadratica con il Metodo delle Celle essendo u la componente del vettore spostamento in un punto P x y z secondo li esimo asse del riferimento globale Adottando il sistema di riferimento locale di coordinate affini di cella introdotto nel paragrafo 3 1 1 1 varr l espressione eo pote I f 3 37 dove il vettore dei coefficienti viene determinato imponendo che la funzione approssimante colga i valori nodali aic C Uic 3 38 essendo aic il vettore dei coefficienti del polinomio interpolante il campo della com ponente u dello spostmento all interno della celle c e C la matrice definita nella 3 20 e uic il vettore delle i esime componenti di spostamento nodale per i nodi ordinati della cella c Con il processo gi illustrato per il campo scalare si perverr pertanto per cia scuna componente di spostamento alle espressioni dine MESE ie eC ag AE He LL LAME PLC 3 39 wlEng 1 CE EC PR CVC we che descrivono completamente la stima del vettore spostamento
227. tta slice In una 11 1 La creazione di modelli numerici subject specific di segmenti ossei Figura 1 2 Schema di funzionamento di uno scanner TC situazione ideale l intensit di un raggio che attraversa una regione a densit co stante si attenua in modo esponenziale lungo il suo percorso Il logaritmo del rapporto fra intensit incidente e intensit trasmessa corrisponde all integrale di linea del coefficiente di attenuazione lineare lungo il percorso all interno del mez zo In un sistema CT ideale la ricostruzione dell immagine ottenibile con una semplice procedura di ricostruzione della proiezione basata su un algoritmo di convoluzione dei segnali di assorbimento misurati Di fatto nei sistemi CT reali sono presenti una serie di imperfezioni che devono essere corrette per minimizzare gli artefatti nella ricostruzione delle immagini correzioni di offset e di guadagno per eliminare gli errori tra i sensori procedu re di correzione dell indurimento del fascio per minimizzare gli effetti dovuti alla dispersione dello spettro di energia del fascio di raggi X altre tecniche di compensazione permettono di eliminare artefatti dovuti alla specifica geometria dello scanner Nell immagine ricostruita le differenze di contrasto sono dovute alle diffe renze del coefficiente di attenuazione lineare nei tessuti attraversati dai raggi X lartefatto che si produce quando presente una brusca variazione nella densit del me
228. un pixel nel piano rispetto al datset CT che costituisce invece il riferimento nativo del modello MCM Questo fatto pu quindi generare un errore di lettura delle informazioni fatto particolarmente influente nelle regioni di confine 169 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle e I due modelli hanno risoluzioni non confrontabili il modello MCM co struito alla risoluzione nativa del dataset laddove il modello FEM carat terizzato da elementi di dimensione variabile tendenzialmente crescente dalla superficie all interno gli elementi del modello FEM mediano infor mazioni materiali HU e quindi E associati ai punti contenuti su un vo lume di dimensioni consistenti laddove il calcolo attuato dall MCM consi dera valori nodali alla risoluzione nativa valori del polo e dei 6 satelliti Nelle regioni di tessuto spongioso il tessuto osseo organizzato in strut ture denominate trabecole caratterizzato da un elevata porosit che an che laddove abbia scala inferiore alla risoluzione del dataset determina co munque valori locali delle propriet meccaniche ridotti con conseguente sovrastima locale delle deformazioni nel modello MCM Va a tale proposito rilevato come al di l del confronto in oggetto la procedura MCM possa paradossalmente essere pi pertinente al vero quando predica localmente nelle regioni a minore densit puntuale maggiori deformazioni il modell
229. vengono di molto chiarite da qualche esempio il gi ci tato caso del flusso ad esempio di calore attraverso una superficie per defi nizione la quantit di calore che fluisce attraverso una superficie assegnata in un dato intervallo di tempo Il flusso dunque una grandezza naturalmente associa ta agli elementi di superficie e agli intervalli di tempo Si tratta ora di stabilire se l orientazione degli elementi di superficie sia interna o esterna si osserva allora che il flusso indica una quantit in transito attraverso al superficie cambiando il verso positivo di attraversamento della superficie il valore del flusso cambia di segno Poich indicare un verso di attraversamento positivo significa precisare un orientazione esterna si ha evidentemente che il flusso riferito ad elementi di superficie orientati esternamente Per proseguire nell ambito della termostatica risulta evidente che la temperatura una variabile associata ai punti e che la differenza di temperatura di conse guenza associata a segmenti cambiando il verso di percorrenza del segmento la differenza di temperatura cambia segno l orientazione dei segmenti di riferimen to quella interna Analoga situazione si ha per le attribuzioni agli elementi temporali Una gran dezza pu essere associata agli istanti oppure ad intervalli di tempo nella prece dente sezione 2 3 si parlato di istanti e di intervalli primali e duali gli elementi primali dotati
230. y fy 4 43 Tg a bzr cy dz exy fy gx y hry ix ove si inteso indicare con il pedice il numero dei nodi satellite Come regione tributaria del nodo 7 stato assunto un cerchio di raggio pari a met del valore medio della distanza media dei satelliti dal nodo polo 66 Le espressioni della i esima riga della matrice fondamentale esprimente il flusso netto attraverso la superficie di confine della regione tributaria del polo di cui alla formula 4 36 sono rispettivamente kpa tkrR 0 0 0 1 1 MG kps tkr R 0 0 0 1 0 1 Mz 4 44 kpg tkrR 0 0 0 2 0 2 0 0 R Me La stima dell ordine di convergenza ossia del decrescere dell errore al diminuire di una misura significativa della distanza dei punti stata effettuata con riferi mento ad una griglia regolare di punti situazione in cui il passo di griglia d una misura esatta della dimensione associata alla distribuzione di punti La tabella 4 4 riporta i valori ottenuti per il parametro 6 con passi della griglia man mano ridotti passo di errore medio griglia 4satelliti 5 satelliti 8 satelliti 0 5 0 00376808 0 00376808 0 00384507 0 2500 0 00078737 0 00078737 0 00062571 0 1250 0 00017645 0 00017645 0 00013119 0 0625 0 00004159 0 00004159 0 00003041 Tabella 4 4 Valore dell errore medio per i punti interni al dominio al variare del passo della regolare di punti con differenti composiz
231. za la procedura viene qui riassunta dedicando maggiore attenzione alle operazioni specifiche del problema ora trattato considerato un generico nodo P interno al dominio assunto come nodo polo e costruire la mesh locale di celle primali con un conveniente numero di nodi satellite scelti tra quelli prossimi circostanti in modo da ottenere un insieme di celle triangolari adiacenti tutte aventi in comune il nodo polo P e aventi adeguata regolarit di forma come suggerito dalla figura 4 3 e definire la regione tributaria del nodo polo provvedendo alla costruzione del complesso duale locale come suggerito dalla figura 4 4 e cella per cella si calcolano le forze agenti sulla superficie di confine tra la cella duale del nodo polo e il resto della cella queste grandezze saranno 89 4 L approccio meshless con il Metodo delle Celle funzione dei valori delle componenti dello spostamento dei vertici della cella esaminata Con riferimento alle notazioni di figura 4 9 la forza agente sull elemento di superficie piana orientata A esprimibile come ie ely a D 48 Cc essendo Figura 4 9 Forza agente su una superficie orientata A nella generica cella c orientazione dei vettori area a puro titolo di esempio stata assunta l esatta geometria della cella 1 di vertici P 5 59 A e A le componenti del vettore area A D la matrice del legame costitutivo che per uno stato piano di tensione in un materiale
232. zionario 4 2 2 Sistemi elastici piani in regime di piccoli spostamenti 4 3 Interpolazione locale con funzioni polinomiali 4 3 1 Campo scalare stazionario o oo 4 3 2 Sistemi piani deformabili 0 4 4 Interpolazione locale con funzioni di base radiale 4 4 1 Interpolazione RBF di campi scalari 4 4 2 Interpolazione RBF dicampivettoriali 5 Applicazione esplorativa della formulazione Meshless a celle locali con il Metodo delle Celle 5 1 La selezione delle metodologia meshless applicata 9 2 IieSEESpioranygo arisa litio heart eg 5 2 1 L esperimento di riferimento sx daria Paga aa 5 2 2 IdatasetCI n asolo io ere ae Mec de ae Rod 5 2 3 Il modello FEM subject specific 5 2 4 Il modello MCM i sia aa 5 2 5 La definizione delle propriet materiali 5 2 6 Determinazione dell accuratezza dei modelli numerici 5 5 Risultati aa ws eh we oe eA e a o ra i 6 Validazione della formulazione Meshless a celle locali con il Metodo delle Celle 6 1 L esperimento di riferimento 6 1 1 Le misure sperimentali x inratiso a erano a ria Fe 6 1 2 6 1 3 6 1 4 6 1 5 emda PEE i lie Bete BS le eS kee dti Acquisizione CT registrazione spaziale e proprieta mecca ANE Pete aae I PIERRE RE PI model MEM elenca dio pda o Definizione dell accuratezza 00 4 6 2 Lecriticitaemerse
233. zioni con 15 estensime tri posizionati in punti anatomicamente significativi Accanto alle misure speri mentali un ulteriore elemento di confronto e verifica stato reso disponibile dai modelli FEM realizzati secondo una procedura validata e replicanti le medesime prove meccaniche Questo studio numerico sperimentale 86 costituisce a tut t oggi il pi consistente benchmark di confronto relativamente allo studio della deformazione nei segmenti di ossa lunghe a livello di organo Le analisi condotte con modelli MCM costruiti sulle immagini cliniche dei femori studiati hanno avuto quindi disponibile come per il test preliminare del capitolo 5 il doppio confronto con le misure sperimentali e con le predizioni dei modelli FEM Un ulteriore elemento di valutazione stato assunto con un confronto locale sistematico con uno dei modelli FEM in una delle condizioni di carico confron tando i valori principali di deformazioni ai centroidi degli elementi FEM con i valori estratti dal modello MCM nei medesimi punti 6 1 L esperimento di riferimento L esperimento stato condotto su quattro paia di femori umani provenienti da donatori maschi adulti deceduti per infarto del miocardio e senza pregresse pa tologie all apparato muscolo scheletrico L indagine non distruttiva DEXA Dual linteso quale implementazione di una procedura cognitiva di previsione e descrizione del comportamento di un sistema 150 6 1 L esperimento di riferiment
234. zioni dei nodi satellite attorno al polo non rendono infatti possibile l interpolazione locale con la funzione poli nomiale adottata Che possano sussistere delle configurazioni geometriche ina deguate ad individuare una funzione interpolante in grado di cogliere i valori ai nodi un fatto intuibile se la funzione interpolante fosse di tipo lineare occorre rebbero evidentemente tre punti non allineati per poterla appoggiare ai valori nodali Talune disposizioni di punti rendono quindi il problema dell interpolazione mal condizionato in termini generali quando un problema malcondizionato pic cole variazioni nei valori di ingresso si traducono in enormi errori nei valori ot tenuti come risultato Nel caso dell interpolazione si ha la situazione in cui gli errori di macchina nella valutazione ad esempio della posizione geometrica dei punti provocano errori enormi nella stima dei valori della funzione interpolata e in ultima analisi nella stima dei valori delle temperature nodali L inadeguatezza della disposizione della costellazione di punti rispetto alla fun zione interpolante si traduce in un valore numerico elevato di una quantit detta numero di condizionamento della matrice M definito come u R M M 4 38 Pi il problema risulta ben condizionato pi il valore di x risulta essere conte nuto A parit di conformazione della costellazione di punti pu poi presentarsi un altro probl
235. zzo caso tipico di quando nell osso presente un materiale metallico quale un chiodo o una protesi 12 1 3 La definizione di un modello Figura 1 3 Immagine CT di una sezione di cranio i differenti livelli di grigio corrispondono a differenti valori del coefficiente di attenuazione lineare del fascio di raggi X ossia a diversi valori di densit dei tessuti attraversati L immagine prodotta quindi la rappresentazione della distribuzione spaziale dei coefficienti di attenuazione lineare Per i tessuti molli il coefficiente di atte nuazione lineare in prima approssimazione proporzionale alla densit pertan to i valori riportati in un immagine CT sono espressi talvolta anche in termini di densit Le immagini CT costituiscono attualmente la migliore fonte di informazione per la costruzione di un modello subject specific di un segmento osseo questo perch il coefficiente di assorbimento del tessuto osseo molto pi elevato di quello dei tessuti circostanti fatto che produce elevati contrasti dell immagine e permette dunque una pi agevole delimitazione del contorno 1 3 1 2 Le Hounfield Units o CT numbers Dopo l operazione di ricostruzione l immagine CT costituita da un set di nume ri floating point utili per effettuare confronti ma non per la loro rappresentazione Siccome la maggior parte dei sistemi di visualizzazione grafica lavora con valori interi dopo la fase di ricostruzione ma prima della visua

sviluppo, valutazione e applicazione di metodi

Contents

Download Pdf Manuals

Related Search

Related Contents