The loads that bones can sustain strongly depend on the composition and microarchitecture of bone tissue. When studying bone failure, experiments and simulations are often applied. Simulations are based on CT images, so that the individual microstructure of a bone can be used as a geometric basis. However, high resolutions in a range of micrometres, which are required to resolve the tiny sub-structures, lead to immense computational requirements. Thus, most simulations nowadays model bone at the micro-scale as a linear-elastic material. As a result, onset of damage and failure of bone is not directly assessable.
The goal of this thesis is the development of a simulation tool which enables materially nonlin- ear simulations of whole bones at the required high resolutions to take the microstructure into account. Precisely, the thesis follows two main objectives: (1) The extension of an existing, highly efficient, μFE solver with respect to nonlinear material behaviour to enable whole bone simulations at sub-trabecular resolution on widely accessible HPC clusters. (2) Based on exam- ples, investigate how the local results can assist in basic research to gain insight into the failure behaviour of bone.
Starting point of the implementation was the highly efficient, linear μFE solver ParOSol. A simple nonlinear material model consisting of a linear-elastic region, a damage-based degradation region including hardening, and a fracture region was chosen. Other sources of nonlinearity were neglected. The material model was applied locally to each element of the structure. ParOSol was extended for the proposed material model using an incremental-iterative solving procedure. The good parallel properties of the original ParOSol were preserved due to the careful implementation of the nonlinear material behaviour. In a performance study using radius segments, a structure containing more than 600 million elements (2.3 billion degrees of freedom) was success- fully simulated using the new simulation framework. The memory requirements were reduced by almost three orders of magnitude compared to results presented in the literature.
Suitable material parameters were identified using two trabecular biopsies. A validation study using axial tests of trabecular biopsies and radius segments proved that the simple modelling approach was sufficient to achieve excellent correlation compared to experiments. All 20 biopsy samples showed rather diffuse damage patterns before more localised failure occurred. A con- siderable amount of damage was present at the apparent yield point (determined via the 0.2 % strain offset criterion), specifically 2 − 6 % of the elements were damaged. A novel elasticity limit εie was proposed. It indicates the onset of nonlinearity using the number of damaged elements obtained from simulation results. Damage modes were defined to classify damaged elements. Systematic local differences between samples under tension and compression were identified. As expected, samples under apparent tension load failed predominantly due to overloading in local tension. Samples under apparent compression showed dominant compression but also increased tension damage, indicating local compression and bending.
The role of stress relocations at overloaded regions for the apparent failure behaviour of radius segments was studied. A comparison between materially linear and nonlinear simulation re- sults was conducted. (1) Based on a linear-elastic material behaviour, the ultimate load was estimated using the Pistoia criterion with two different approaches: (a) The standard Pistoia criterion where the ultimate load was obtained by scaling the apparent force such that 2 % of the elements sustained more than 0.7 % strain. (b) In a second approach, optimal Pistoia parameters were identified: The load obtained from the linear simulations was fitted to the experimental failure loads, and the number of overstrained elements (strain larger than 0.7 %) was determined. (2) Materially nonlinear simulations were conducted using the previously identified material parameters without any further calibration.
(1a) The standard Pistoia criterion yielded a good correlation, but the ultimate load was severely underestimated. (1b) The optimal Pistoia criterion resulted in much larger overstrained volumes ((39.49 ± 19.45) %) than the standard Pistoia criterion. The optimal Pistoia criterion led to similar overstrained regions compared to the damage regions of the nonlinear simulations. (2) The ultimate force was directly obtained from the nonlinear simulations and showed good agreement with experimental data.
The main progress of this thesis is the extension of an efficient solver to nonlinear material behaviour. A framework was provided which enables solving of huge geometries on HPC clusters, and medium sized geometries in reasonable times on conventional shared memory machines. Another unique aspect of the developed solver is that a failure mechanism and tissue level softening is included. Thus, the maximum point in the apparent force-displacement curve can be directly computed. Additionally, highly damaged and fractured regions can be identified. This makes it possible to investigate the impact of stress relocations and can lead to further insight into the failure behaviour of bone structures.
Die Festigkeit von Knochen hängt stark von ihrer Zusammensetzung und Mikrostruktur ab. Um das Versagen von Knochen zu studieren werden häufig Experimente und Simulationen verwendet. Simulationen basieren meist auf hoch aufgelösten CT-Scans, sodass die individuelle Mikrostruk- tur eines Knochens direkt abgebildet werden kann. Hohe Auflösungen im Mikrometer-Bereich, welche notwendig sind um die kleinen Substrukuren aufzulösen, führen zu einem hohen Rechenbe- darf. Daher wird Knochen auf der Mikroebene meist mit einem linear-elastischen Materialver- halten modelliert. Dies führt jedoch dazu, dass das nichtlineare Bruchverhalten nicht direkt untersucht werden kann.
Das Ziel dieser Arbeit ist die Entwicklung eines Simulationswerkzeugs, welches Finite-Elemente- Simulationen (FE) von ganzen Knochen mit der notwendigen Auflösung und einem nichtlinearen Materialmodell ermöglicht. Diese Aufgabe wurde in zwei Schritten behandelt: (1) Zuerst wurde ein existierendes Programm zur Durchführung von linearen μFE-Simulationen für ein einfaches nichtlineares Materialverhalten erweitert, sodass nichtlineare Simulationen ganzer Knochen auf leicht zugänglichen HPC-Servern möglich sind. (2) Anhand von Beispielen wurde untersucht wie die lokalen Ergebnisse helfen können in der Grundlagenforschung mehr Einsicht in das Bruchver- halten von Knochen zu erhalten.
Als Basis für die Entwicklung wurde der effiziente μFE-Code ParOSol ausgewählt, welcher bisher nur linear-elastisches Materialverhalten erlaubte. Ein einfaches nichtlineares Materialmodell wurde vorgeschlagen, welches basierend auf einer skalaren Schädigungsvariable einen linear- elastischen Bereich, einen Schädigungs- und einen Bruchbereich unterscheidet. Andere Arten der Nichtlinearität wurden vernachlässigt. Zur Implementierung wurde das Materialmodell lokal, d.h. auf jedes Element der Struktur, angewandt. Eine Gleichgewichtslösung wurde über einen inkrementell-iterativen Lösungsprozess ermittelt.
Aufgrund der angepassten Implementierung des neuen Materialmodells in ParOSol konnten die guten numerischen Eigenschaften des ursprünglichen Lösers bewahrt werden. In einer Skalie- rungsstudie wurden erfolgreich Strukturen mit bis zu 600 Millionen Elementen (2,3 Milliarden Freiheitsgraden) simuliert. Die benötigte Menge an Arbeitsspeicher konnte im Vergleich zur Literatur um fast drei Größenordnungen reduziert werden.
Materialparameter für das gewählte Materialmodell wurden an Hand von zwei trabekulären Biop- sien identifiziert. Eine Validierungsstudie an axialen Tests von trabekulären Biopsien und Radius- Segmenten zeigte, dass der einfache Modellierungsansatz ausreicht um exzellente Korrelationen mit Experimenten zu erreichen. Alle 20 Biopsien zeigten eine diffuse Schädigung der Knochen- struktur bevor stark lokalisierte Brüche auftraten. An der apparenten Fließgrenze (ermittelt durch das 0,2 %-Dehnungskriterium) waren schon 2-6 % der Elemente beschädigt. Daher wurde eine alternative Definition zur Ermittlung der Elastizitätsgrenze vorgeschlagen. Beschädigte Elemente wurden in Belastungsmoden eingeteilt. Diese zeigten systematische Unterschiede zwischen Proben unter apparenter Zug- oder Druckbelastung. Wie erwartet, beschädigten Proben unter Zugbelastung hauptsächlich durch lokale Zug-Überbelastungen. In Proben unter Druckbelastung war der Druck-Modus dominant, wobei auch größere Mengen lokaler Zugbelastung auftraten.
Die Rolle von Spannungsumlagerungen in beschädigten Bereichen für das Bruchverhalten von Knochen wurde untersucht. Dazu wurde ein Vergleich zwischen Simulationen mit linearem und nichtlinearem Materialverhalten durchgeführt. (1) Basierend auf linear-elastischem Materialver- halten wurde die Maximalkraft mittels des Pistoia-Kriteriums bestimmt. Zwei Varianten wurden verwendet: (a) Das übliche Pistoia-Kriterium, wo die Maximalkraft skaliert wurde, sodass min- destens 2 % der Elemente über 0,7 % gedehnt wurden. (b) In einer zweiten Variante wurden “optimale” Pistoia-Parameter bestimmt. Dazu wurde die Maximalkraft an die experimentell er- mittelte angepasst und der dazu notwendige Anteil an überdehnten Elementen ermittelt. (2) Das nichtlineares Materialmodell mit den zuvor identifizierten Materialparametern wurde verwendet.
(1a) Das übliche Pistoia-Kriterium führte zu guten Korrelationen, aber die Maximalkraft wurde stark unterschätzt. (1b) Das optimale Pistoia-Kriterium führte zu wesentlich größeren Anteilen an überdehnten Elementen ((39.49 ± 19.45) %). Die Bereiche der überdehnten Elemente des op- timalen Pistoia-Kriteriums waren vergleichbar in Lage und Größe mit den Bereichen beschädigter Elemente der nichtlinearen Simulationen. (2) Die Maximalkraft war direkt aus den nichtlinearen Simulationen ermittelbar und zeigte gute Übereinstimmung mit experimentellen Ergebnissen.
Der wichtigste Fortschritt dieser Arbeit liegt darin, dass ein effizientes Programm um nichtlin- eares Materialverhalten erweitert wurde. Einerseits ist es möglich, große Geometrien auf gängi- gen HPC-Servern zu simulieren. Anderseits können mittelgroße Geometrien auf konventionellen Shared-Memory-Computern in akzeptablen Rechenzeiten gelöst werden. Ein weiterer besonderer Aspekt des Programms ist, dass ein Bruch- und Beschädigungsmechanismus inkludiert wurde. Dadurch kann die maximale Traglast direkt aus einer Simulation bestimmt werden. Es lassen sich sowohl fortschreitende lokale Schädigung als auch der makroskopische Bruch von Strukturen verfolgen. Außerdem sind tiefere Einblicke in das Schädigungsverhalten von Knochen möglich.