An extended biphasic description of the inhomogeneous and anisotropic intervertebral disc

Abstract

It is the aim of this contribution to develop a finite element model, which is as simple as possible, but at the same time complex enough to capture many of the occurring tissue properties of the intervertebral disc (IVD). In order to better understand these properties from an engineering point of view, the needed basic anatomical knowledge is briefly reviewed in the beginning of this treatise, thereby addressing the lumbar spine with focus on the IVD and its material properties. In particular, the IVD appears as the largest avascular part of the body and its microstructure leads to an electro-chemically active material with anisotropic, inhomogeneous and strongly dissipative behaviour.

In the following main part of this work, the complete continuum-mechanical modelling process is extensively discussed as well as the numerical treatment of the resulting governing equations. Starting from the thermodynamically consistent Theory of Porous Media (TPM), two phases and three components are introduced for the description of IVD tissue. In particular, this is the extracellular matrix (solid skeleton) carrying fixed negative charges which is saturated by a pore fluid consisting of a solvent (liquid) as well as anions and cations of a dissolved salt. Following the idea of superimposed continua, an individual motion function is introduced for each of the constituents, whereas the components of the pore fluid are always expressed relative to the deforming solid skeleton. In order to capture the finite kinematics of the inelastic solid skeleton, its deformation gradient is multiplicatively split into inelastic and elastic parts. Next, the materially independent balance equations are derived from the respective master balances and accustomed to the soft biological tissue under study. In order to keep the resulting set of equations as simple as possible, while still keeping the ability to reproduce osmotic effects, an assumption according to Lanir is made. In this context, the tissue is regarded to be always immediately in electro-chemical equilibrium, which allows to describe the electro-chemically active tissue using only an extended biphasic model. Applying van't Hoff's law finally allows to compute the occurring osmotic pressure as a function of the solid displacement.

Moreover, in order to characterise the inhomogeneous anisotropic and viscoelastic solid skeleton as well as the viscous pore fluid, several constitutive equations need to be formulated, thereby depending on a thermodynamically admissible set of process variables. Herein, the endangerment of postulating nonphysical constitutive assumptions is avoided by strictly following the restrictions resulting from the evaluation of the entropy inequality. Finally, the chosen constitutive functions of the solid skeleton are based on Ogden-type strain energy functions, which automatically include several simpler material laws. The viscoelastic contribution is based on a generalised Maxwell model which is dominated by the concept of internal variables with linear evolution equations. Finally, the superimposed dissipative effect of the viscous pore fluid is captured using the famous Darcy filter law.

As a last step, the applicability of the derived model is proven with realistic computations of the IVD. Herein, the resulting set of governing partial differential equations is discretised in time and space using the finite difference method and the mixed finite element method, respectively. The theoretically introduced material parameters are determined using experimental data as well as material parameters obtained from a vast collection of related literature sources. Since many parameters appear in a coupled manner, their identification is often only possible via inverse computations. Following this, a numerical sensitivity analysis is carried out yielding an indication for the relevant parameters in experiments concerning a motion segment in a short-duration compression-flexion experiment as well as in long-term loading situations. Subsequently, the efficiency of the implementation is demonstrated by a parallel simulation of a lumbar spine segment carried out on 84 processors simultaneously, thereby exhibiting almost one million degrees of freedom.

Das Ziel der vorliegenden Arbeit ist die Entwicklung eines Finite-Elemente-Modells der menschlichen Bandscheibe. Dieses soll auf der einen Seite so einfach wie möglich gehalten werden, auf der anderen Seite aber komplex genug sein, um die relevanten Gewebeeigenschaften der Bandscheibe abzubilden. Um diese Eigenschaften als Ingenieur besser verstehen zu können, werden die anatomischen Grundkenntnisse der Lendenwirbelsäule einleitend bereitgestellt, um anschließend auf die Zusammensetzung der Bandscheiben detailliert einzugehen. Dabei entpuppt sich die Bandscheibe als größter Vertreter avaskulärer Gewebe im Körper und repräsentiert ein elektro-chemisch aktives, anisotropes, inhomogenes und stark dissipatives Material.

Der nun folgende Kern der Arbeit befasst sich mit dem vollständigen Spektrum der kontinuumsmechanischen Modellbildung sowie deren numerische Umsetzung. Ausgehend von der thermodynamisch konsistenten Theorie Poröser Medien (TPM) werden zwei Phasen und drei Komponenten für die Beschreibung des Bandscheibegewebes eingeführt. Im Einzelnen sind dies der Festkörper, bestehend aus einer extrazellulären Matrix (Solid) mit anhaftenden gebundenen Ladungen, und ein Porenfluid bestehend aus einem Lösungsmittel (Liquid) und darin gelösten Kationen und Anionen. Im Rahmen superponierter Kontinua werden individuelle Bewegungsfunktionen für die Konstituierenden eingeführt, wobei die Fluidkomponenten immer relativ zum sich deformierenden Festkörperskelett ausgedrückt werden. Um die finite Kinematik des inelastischen Festkörperskeletts zu beschreiben wird dessen Deformationsgradient multiplikativ zerlegt. In einem nächsten Schritt werden die zur Modellierung benötigten Bilanzgleichungen anhand einer Masterbilanz hergeleitet und auf das zu beschreibende Gewebe angewendet. Damit das Endergebnis relativ simpel bleibt, jedoch gerade noch osmotische Effekte abbilden kann, wird eine Annahme nach Lanir getroffen, die von instantanem elektro-chemischen Gleichgewicht im biologischen Gewebe ausgeht. Diese Vereinfachung erlaubt es wiederum ein elektro-chemisch aktives Material mit einem erweiterten Zweiphasenmodell zu beschreiben, da mit Hilfe des van't Hoff'schen Gesetz der osmotischen Druck im Gewebe rein als Funktion der Festkörperverschiebung berechnet werden kann.

Damit das viskose Porenfluid sowie das inhomogene, anisotrope und viskoelastische Festkörperskelett beschrieben werden können, bedarf es einer konstitutiven Grundauswahl an Antwortfunktionen, die von thermodynamisch zulässigen Prozessvariablen abhängen müssen. Um bei der Konkretisierung dieser Antwortfunktionen keine unphysikalischen Annahmen zu treffen, werden die maßgeblichen Restriktionen aus der Entropieungleichung herausgearbeitet. Die gewählten Konstitutivgleichungen des Festkörperskeletts basieren auf Verzerrungsenergiefunktionen vom Ogden-Typ, die somit automatisch auch einfachere Gesetze beinhalten. Der viskoelastische Anteil bezieht sich auf ein generalisiertes Maxwell Modell und wird durch das Konzept der internen Variablen mit linearen Evolutionsgleichungen geprägt. Darüber hinaus liefert das viskose Porenfluid einen überlagerten dissipativen Effekt, der durch das Darcy'sche Filtergesetz berücksichtigt wird.

Um die Anwendbarkeit des Modells zu demonstrieren werden die resultierenden partiellen Differentialgleichungen im Rahmen der gemischten Finite-Elemente-Methode numerisch diskretisiert und mittels geeigneter dreidimensionaler Anfangs-Randwert-Probleme einer Anwendung zugänglich gemacht. Die benötigten Materialparameter werden anhand von Versuchsergebnissen bestimmt, die der entsprechenden Literatur entnommen wurden. Des Weiteren wird eine Sensitivitätsanalyse der involvierten Parameter durchgeführt, um den Einfluss einzelner Parameter auf das gesamte Deformationsverhalten eines Bewegungssegments während kurzzeitiger Flexion und langzeitiger Kompression zu bestimmen. Um die Effizienz des Modells abschließend zu verdeutlichen, wird mit einer parallelen Lösungsstrategie ein großskaliges dreidimensionales Problem der Lendenwirbelsäule mit nahezu einer Million Freiheitsgrade auf 84 Prozessoren parallel gelöst.