DE102005035181A1

DE102005035181A1 - Verfahren zur Darstellung von dreidimensionalen Strömungsdaten und Wandbelastungen in fluidführenden Gefäßen, insbesondere Blutgefäßen, unter Verwendung der Lattice-Boltzmann-Methode

Info

Publication number: DE102005035181A1
Application number: DE102005035181A
Authority: DE
Inventors: Kamen Dr. Beronov
Original assignee: Friedrich Alexander Universitaet Erlangen Nuernberg
Current assignee: Friedrich Alexander Universitaet Erlangen Nuernberg
Priority date: 2004-07-27
Filing date: 2005-07-27
Publication date: 2006-03-23
Anticipated expiration: 2025-07-28
Also published as: DE102005035181B4

Abstract

Dreidimensionale (3-D) Grauwertdaten aus MRT oder CT (Röntgen-Tomographie) werden mit Varianten geeigneter Algorithmen für die medizinische und allgemeine Bildgebung binär segmentiert. Entstandene 3-D-Gefäßgeometrien werden für eine aufgelöste numerische Simulation ihrer Durchströmung auf kubischen Rechengittern mit Lattice-Boltzmann-Strömungslösern verwendet. Benötigte Randbedingungen werden aus Datenbanken oder aus individuellen nichtinvasiven Messungen extrahiert. Wenn die Wanddynamik wichtig ist, wird sie mitsimuliert. Im Patienten vorhandene Implantate, Vasoplastik, vorhandene Plaquen etc. werden identifiziert und in die Simulationsgeometrie eingebracht. Gleichfalls werden auch virtuelle (im Patienten noch nicht entstandene) Geometrieänderungen, z. B. einzubringende Implantate, zu entfernende Gefäßteile etc., an die zu berechnende Geometriebeschreibung angebracht, um bei der medizinischen Planung von entsprechenden operativen Eingriffen eine vorzeitige präzise Einschätzung dessen Auswirkungen zu ermöglichen. Implantate, dessen Feinstruktur auf dem Rechengitter der Strömungssimulation nicht unaufgelöst werden kann, werden als elastische, poröse Körper modelliert.

Description

Die Erfindung betrifft ein Verfahren zur Darstellung von zeitabhängigen, räumlich dreidimensionalen Strömungsdaten in fluidführenden Gefäßen unter Verwendung der Lattice-Boltzmann-Methode, bei dem dreidimensionale Gewebestruktur-Daten auf regelmäßigen Datengittern (Geometrie-Voxeldaten) des darzustellenden Gefäßbereichs einschließlich dessen Wandung, d.h. deren geometrische Struktur darstellende Dichtedaten, verwendet werden und ebenfalls regelmäßige isotrope Rechengitter (Voxel-Gitter) für die Gefäßdarstellung und für die Bestimmung des Strömungsverlaufs ausgewählt werden, wonach der Verlauf der Fluidströmung in zeitlichen Schritten bestimmt wird.
Die Ermittlung eines Strömungsfeldes, insbesondere in einem Gefäßsystem, ist sehr aufwendig. Für die hydrodynamische Bestimmung müssen die Geometrie des durchströmten Gebiets, die Ein- und Ausström-Randbedingungen und die Bedingungen an dessen Wänden gegeben sein. Die Ein- und Ausström-Randbedingungen umfassen zeitabhängige Angaben zum Druck und, insofern bekannt, zu den Geschwindigkeiten und im Falle turbulenter Strömungen auch noch zu den Spannungen. Diese Angaben müssen die räumliche und zeitliche Abhängigkeit der Strömung in einer ausreichenden Genauigkeit enthalten.
Blutströmungen durch Adern, insbesondere Aortenzweige und weitere größere Arterien, sind von einer starken pulsierenden Kom ponente und von komplizierten Strömungsmustern in der Nähe von Bifurkationen, Krümmungen sowie Malformationen (Erweiterungen oder auch Verengungen) geprägt. Die Gefäßwände sind schnellen elastischen Deformationen im Herzpuls-Rhythmus ausgesetzt, was für die Hämodynamik je nach Lage und Zustand der Gefäße eine wesentliche Rolle spielen kann. Die Modellierung der mechanischen Eigenschaften dieser Gefäßwände ist nicht eindeutig und es ist besonders schwierig, individuelle und ortsspezifische Korrekturen einzuführen, da die gekoppelte Bestimmung der Wand- und Blutdynamik wesentlich aufwendiger als die Behandlung jeder der beiden dynamischen Prozesse allein ist und weil die Beschaffenheit von Gefäßwänden ein allgemein schwieriges, noch offenes Forschungsgebiet darstellt. Die Verifikation durch direkte Messungen ist auch durch die räumlich und zeitlich sehr hohe erforderliche Auflösung erschwert.
Im Bereich der Computational Fluid Dynamics (CFD) verfügbare Standardverfahren zur direkten numerischen Simulation (DNS) inkompressibler Strömungen, meistens unter Anwendung von Finite-Volumen-Verfahren (FVM), ebenso CSM-Verfahren (computational structure mechanics) unter Anwendung von Finite-Elemente-Verfahren (FEM), bei denen eine Diskretisierung mit unstrukturierten Rechengittern des Strömungsvolumens und der Gefäßwände durchgeführt wird, ermöglichen die Durchführung aufgelöster 3D-Simulationen der Durchströmung von größeren Blutgefäßen unter Einsatz insbesondere von CT-Daten. Die Ergebnisse solcher Simulationen sind detaillierte dynamische Informationen über Strömungsstörungen und Wandbelastung, die in den medizinischen Meßdaten selbst nicht enthalten sind. Die Durchführung derartiger Berechnungen mit den benannten Verfahren stellt hohe Anforderungen an die untersuchenden Personen und ist außerordentlich zeitaufwendig, sie kann z.B. für eine Geometrie eine Woche dauern. Eine Vereinfachung der Simulationsverfahren bleibt durch die unentbehrlichen weitgehenden Modellannahmen immer auf spezifische Gefäße, einfache Geometrie der Malformationen und enge Parameterbereiche stark eingeschränkt.
Mit wesentlich aufwendigeren Verfahren, die auf dem Doppler- oder dem Time-of-Flight-Prinzip aufgebaut sind und dabei auf dem allgemeinen mathematischen Apparat der Tomographie basiert sein dürfen, wie z.B. Phasenkontrast-MRI oder z.B. Doppler-Ultraschall-Anemometrie, können auch einzelne Geschwindigkeitskomponenten und Durchschnittsgeschwindigkeiten, also ein Teil der relevanten kinetischen Information, aufgenommen werden. Die Auflösung ist in meisten Fällen gering, so daß nur das Gesamtvolumen, das ein bestimmtes Lumen pro Zeiteinheit durchströmt, mit hinreichender Zuverlässigkeit gemessen werden kann. Lokale Störungen, starke Gradienten und Turbulenz können nicht quantitativ zuverlässig identifiziert werden, obwohl ihre Spuren und Einflüsse mit manchen Verfahren entdeckt werden können. Eine allerdings eingeschränkte 3D-Information mit guter Auflösung liefert auch das nichtinvasive DSA-Verfahren (Digital Substraction Angiography), ein Standardverfahren für die Ermittlung der kinetischen Blutstrominformation vor und während der Gefäßchirurgie. Dagegen kann dynamische Information, d.h. die zeitabhängige Kraftbilanz, bei Blutgefäßen entweder nur indirekt oder wie bei Blutdruck auch direkt aber dann invasiv und nur für einen Teil der Kräfte gemessen werden. Kein zur Zeit verfügbares Meßverfahren, weder indirekt noch in situ, ist in der Lage, 3D-aufgelöste dynamische Messungen zu liefern.
Bei Blutgefäßen liegt das medizinische Hauptinteresse nicht so sehr auf der Geschwindigkeitsverteilung, als auf der damit verbundenen Verteilung von Druck und Scherspannung an der Gefäßwänden. Es ist der Medizin wohlbekannt geworden, daß große Abweichungen der Scherspannungen zur Restrukturierung, Remodellierung und möglicherweise zur Beschädigung der Gefäßwände, beispielsweise im Zusammenhang mit Sklerose, Stenose, hohem Blutdruck und ähnlichen Erkrankungen, führen, die das Leben der Betroffenen behindern und auch direkt bedrohen. Für die Diagnostik dieser Art von Gesundheitsbedrohungen kann aber auch die Analyse der Geschwindigkeitsverteilungen selbst gut gebraucht werden. Noch wichtiger ist es, daß die Rolle des Drucks (z.B. bei Aneurysmen offensichtlich ein Hauptfaktor, der zu Rupturen führen kann) und des Wirbelfeldes wenig bekannt ist und deswegen leichter und besser zu analysieren sein soll als bisher.
Zunehmend wird in der Fluidmechanik die Bestimmung komplexer Strömungen unter Einsatz der Lattice-Boltzmann-Methode (LBM) durchgeführt, die auf einer regelmäßigen Gitterstruktur basiert. Die Zellen eines solchen Gitters können mittels Hochleistungsrechnern vergleichsweise schnell verarbeitet werden. So wurde beispielsweise von T. Zeiser in "Numerische Simulation reaktiver Strömungen in komplexen Geometrien mit der Lattice Boltzmann Methode", LSTM-Jahresbericht, 2001 eine Geometrie-Modellierung vorgestellt, die auf einer Abbildung realer Proben aus Rohrbündelreaktoren mittels 3D-Röntgen-CT auf ein kartesisches Voxelgitter beruht.

In der DE 100 50 063 A1 ist ein Verfahren der eingangs genannten Art zum rechnergestützen Ermitteln eines Strömungsfeldes in dem Atemweg eines Tieres vorgeschlagen worden, um die Strömungseigenschaften der ein- und ausgeatmeten Luft durch den Nasenraum eines Patienten zu ermitteln und bei einer strukturellen Änderung des Atemwegs die voraussichtlich entstehenden Strömungseigenschaften innerhalb des veränderten Atemwegs zu ermitteln. Dabei werden die Wände des Nasenraums durch Segmentierungstechniken zu zusammenhängenden CT-Bild-Elementen bis zur Genauigkeit der Voxel-Geometrie erkannt. Für jedes CT-Bild wird ein kartesisches Gitter gebildet und so das Bild in einzelne Gitterelemente unterteilt, wobei unter diesen Gitterelemente unterschieden werden, die sich innerhalb des durch die Struktur definierten Innenraums, also im Atemweg, befinden und für diese Gitterelemente ein Strömungsfeld ermittelt. Aus der gegebenen Beschreibung des Verfahrens ist zu schließen, daß nur eine einfache, Voxel-basierte Behandlung der Haftrandbedingungen an den Wänden der Atemwege, z.B. mit dem Bounce-back-Verfahren, in Frage kommt und die Krümmung der Wände nur stufenweise, der Auflösung der CT-Daten entsprechend dargestellt werden kann. Besonders bei höheren Strömungsgeschwindigkeiten wird dies zu einer großen Fehlerquelle. In der medizinischen Praxis ist das gesamte Verfahren zur Simulation des Atemvorgangs jedenfalls derzeit nicht anwendbar, da die zur zuverlässigen Berechnung der verhältnismäßig schnellen Luftströmung (hohe Reynoldszahlen) benötigten Rechenzeiten immer noch sehr groß sind.

Zur Behandlung zerebraler Aneurysmen mit einem porösen Stent ist von B. Chopard, "A Lattice Boltzmann Study of Blood Flow in Stented Aneurysmen", SCS Seminar – Computational Hämodynamics, Universität Amsterdam, 13. Oktober 2003 als Maßnahme für die Strömungsreduktion eine Darstellung der Blutströmung mittels der Lattice-Boltzmann-Methode vorgeschlagen worden. Dabei handelte es sich aber um eine zwei-dimensionale Simulation, die für keine reale klinische Anwendung eine genügende Erfassung der individuell vorliegenden Gefäßgeometrie ermöglicht und sogar grundsätzliche Strömungseffekte, die mit dem dreidimensionalen Charakter der Blutströmung zusammenhängen, falsch einschätzt. Bisher sind bei Anwendungen der Lattice-Boltzmann-Methode auf Blutbahnen jedoch nur ein einziges quaderförmiges Gitter und ein einfaches Rechenschema für die Simulation der Strömungsvorgänge eingesetzt worden.

Die US 2003/0060988 A1 befaßt sich mit der Analyse des Einfüllens in Gießformen, wobei die Lattice-Boltzmann-Methode aufgrund der schnellen Rechenzeit und ihrer Eignung für eine große Anzahl von Diskretisationspunkten angewendet wird. Das darin beschriebene Rechenverfahren ist allerdings nur für Strömungen mit relativ niedrigen Geschwindigkeiten gut geeignet, z.B. aufgrund der geringeren Anzahl von Freiheitsgraden und der Herleitung der Randbedingungen, es ist sehr kompliziert und es ist an sich nicht geeignet, fortgeschrittene Blutrheologie-Modelle zu behandeln.

Besonders die Herz-, Gefäß- und Neurochirurgie macht von neueren technischen Entwicklungen und immer genaueren, effizienteren vorzugsweise minimal-invasiven Verfahren Gebrauch, um dreidimensionale (3D) medizinische Daten auszuwerten und zu interpretieren. Die zur Zeit für die klinische Praxis üblichen Tomo graphiedaten, die üblicherweise auf Ultraschall-, Magnetfeld- und Röntgenstrahl-Messungen basieren (Die Erfindung soll aber nicht auf diese beschränkt sein.), enthalten Dichte-Informationen, die danach durch sogenannte Segmentierung vorwiegend in geometrische Information umgewandelt wird. Die Auflösung dieser Verfahren verbessert sich in der gegebenen Reihenfolge, da die minimal auflösbare Länge proportional zur Strahlungswellenlänge ist. Zur Zeit ist eine isotrope Auflösung von ca. 0,4 mm mit Mehrschicht-Röntgen-CT erreichbar, eine etwas gröbere und dennoch vergleichbare (ca. 1 mm) Auflösung mit MRI. Dies bedeutet, daß Gefäße mit Lumendurchmessern von mindestens 1–2 bzw. 3–4 mm detektiert werden können. Das Verhältnis zwischen der Auflösung, die für die Erkennung einer Ader ausreicht, und der, die für eine Bestimmung der Strömungsdynamik und besonders der Wandbelastung (durch Messungen oder Simulationen) mindestens gebraucht wird, ist etwa 1:10.

Es ist also noch nicht möglich, detaillierte Analysen des Blutstromverlaufs und von dessen Zusammenhang mit der Entwicklung von Gefäßanomalien nur anhand medizinischer Meßverfahren durchzuführen. Andererseits hat eine Vielfalt von medizinischen Untersuchungen verschiedene Aspekte festgestellt, bei denen Änderungen im Blutstromverlauf die Entstehung und Weiterentwicklung von Gefäßmalformationen provozieren. Zum Verständnis und zur Behandlung von Gefäßbeschädigungen sind detaillierte Kenntnisse über den Blutströmungsverlauf und des Einflusses mechanischer Belastungen auf die Deformation wichtiger Blutgefäße notwendig. Wünschenswert wäre daher ein Verfahren, das gleichzeitig effizient und einfach Strömungsverhältnisse nachbilden und detaillierte hämodynamische Information erstellen kann.

Um individuell spezifische Merkmale der Gefäßgeometrie erfolgreich einzubeziehen, ist von A. Taylor, M. T Draney, J. P. Ku, D. Parker, B. N. Steele, K. Wang, and C. K. Zarins, "Predictive medicine: computational techniques in therapeutic decisionmaking", Computer Aided Surgery 4, 231 bis 247, 1999 eine Methode vorgeschlagen worden, mit der der Querdurchmesser des Lumens größerer Arterien anhand von MRI-Tomographie-Daten automatisch bestimmt und diese Information für 1D-Modellberechnungen des Blutstroms (nur entlang der Gefäßachse) durch diese Arterien benutzt wird, auch wenn Implantate (incl. Bypässe) vorhanden sind. Bisher ist diese Methode bei der AA (ascending aorta) in Tieren und Menschen getestet worden, also in einem Gefäß, das ausnahmsweise gerade und lang ist. Vergleiche mit MRI-basierter Anemometrie, veröffentlicht in B. N. Steele, J. Wan, J. P. Ku, T. J. R. Hughes und C. A. Taylor, "In vivo validation of a one-dimensional finite-element method for predicting blood flow in cardiovascular bypass grafts", IEEE Trans. Biomed. Eng. 50(6), 649 bis 655, 2003 und mit detaillierten numerischen 3D-Simulationen, veröffentlicht in V. Favier and C. A. Taylor, "Image-based modeling of blood flow in a procine aorta bypass graft", Ann. Res. Briefs, Center for Turbulence Research, Stanford Univ., 2002 haben gezeigt, daß manche gemittelte Größen gut reproduziert werden können, daß die räumliche Verteilung der mechanischen Spannungen und der Blutgeschwindigkeiten jedoch nicht vorhergesagt werden kann.

Die Simulation der Blutströmung mit Berücksichtigung realer Gefäßgeometrien, einschließlich beispielsweise von Implantaten wie Stents oder Aneurysmen, und deren Vernetzung (siehe z.B. M. S. Juchems, D. Pless, T. R. Fleiter, A. Gabelmann, F. Liewald, H. J. Brambs, A. J. Aschpoff, "Blutflußsimulationen mittels Computational-Fluid-Dynamics an aus CT-Daten rekonstruierten Aorten-aneurysmata vor und nach Stent-Graft Implantation", Fortschr. Röntgenstr. 176, 56 bis 61, 2004) macht Gebrauch von CT-Daten unter Einsatz von mindestens drei verschiedenen Auswertemethoden incl. Visualisierung. Das deklarierte Endziel der Entwicklung eines Verfahrens für simulationsgestützte medizinische Planung ist auf diesem Wege kaum erreichbar.

Die Diskretisierung des durchströmten Volumens bei dreidimensionaler numerischer Simulation (3D DNS) für CFD-Anwendungen findet entweder auf strukturierten oder auf unstrukturierten Gittern statt. Die Auswahl der Gitterart hat einen großen Ein fluß auf die Komplexität, die Genauigkeit und die Effizienz der Simulationen. Sie bedeutet unter anderem auch eine Auswahl konsistenter diskreter Darstellungen für das Volumen einerseits und für dessen Rand andererseits sowie die Kopplung dieser Darstellungen. Bei Blutgefäßen ist die Geometrie relativ kompliziert und das Verhältnis (in dynamisch relevanten Meßeinheiten) von Wandfläche zu Volumen ist relativ hoch, so daß die Notwendigkeit effizienter und zuverlässiger Kopplungsverfahren ausgeprägt ist.

Stukturierte Gitter, auf denen beispielsweise die Lattice-Boltzmann-Methode basiert, entstehen nach einer (groben) a priori Verteilung des Gesamtvolumens in benachbarte oder auch teilweise überlappende Einzelgebiete, die jeweils einem einfachen Körper, z.B. einem kubischen Würfel, diffeomorphisch sind. Dieser einfache Körper wird mit einem regelmäßigem, meistens orthogonalen Gitter (3D räumlich) bedeckt. Für solche Gitter stehen hochgenaue und hocheffiziente Algorithmen zur Verfügung, wie z.B. FFT, kompakte Finite-Differenzen, sowie polynomiale Zerlegungen hoher Ordnung, die z.B. bei der Konstruktion von p-FEM oder Spektral-Elemente-Verfahren verwendet werden. Dafür ist aber der Aufbau von strukturierten Rechengittern für komplexe Geometrien sehr aufwendig oder gar unmöglich ohne vereinfachende Annahmen. Auf unstrukturierten Gittern ist dagegen der Aufbau des Gitters weitgehend automatisierbar, da die restriktive Verteilung des 3D-Raumes vermieden werden kann. Dadurch können auch sehr komplexe Geometrien behandelt werden. Dies führt allerdings im Endeffekt zu höheren Kosten (zum Teil wegen niedriger Genauigkeit, wegen höherer algorithmischer Komplexität, usw.) während der Simulation.

Die Identifikation des Strömungsgebiets bei Simulationen der Hämodynamik in Blutgefäßen basiert auf einer Identifikation der Gefäßwände aus insbesondere CT-Daten, die Grauwert-Daten auf einem einfachen strukturierten Gitter darstellen und dessen Bearbeitung auf der Annahme ihrer ausreichenden Auflösung und in meisten Fällen auch ihrer räumlichen Glattheit beruht. Stan dardverfahren der Bildgebung wie z.B. der Marching-Cubes-Algorithmus können sehr schnell aus diesen Daten eine diskretisierte, approximierende Darstellung der Wandgeometrie, etwa als unstrukturierte 2D-Gitter, erstellen. Die Qualität diese Darstellung ist notwendigerweise eine Funktion der CT-Daten-Qualität, so z.B. wird bei einem Bias in der CT-Aufnahme ein lokaladaptives Schwellwert-Verfahren wie aus R. C. Gonzalez and R. E. Woods, "Digital image processing", Prentice Hall, 2002 notwendig. Die entstehende 2D-Triangulation der Wandoberflächen kann bei DNS auf unstrukturierten Gittern schon direkt für eine automatische 3D-Gittergenerierung benutzt werden. Bei strukturierten 3D-Gittern ist die Triangulation ein allein stehendes Objekt, dessen Verhältnis zum 3D-Rechengitter durch Kopplungsschemen verschiedener Art, Genauigkeit, Komplexität und Effizienz gestaltet werden kann. Unstrukturierte Gitter sind an die Randgeometrie angepaßt und müssen sich daher bei einer Bewegung der Wände mitverändern, was zu aufwendigen globalen 3D-Rechnungen, zur Notwendigkeit der Kontrolle des Deformationsgrades, möglicherweise zur Erzeugung und Löschung von Gitterelementen sowie zur Entstehung von mehreren Fehlerquellen, die mit diesen Operation verbunden sind, führt.

Es besteht aber auch eine praktische Möglichkeit, das 2D-Gitter mit einem Lagrange'schen Verfahren bei jedem Zeitschritt zu verändern, ohne das 3D-Rechengitter modifizieren zu müssen. Es ist bemerkenswert, daß eine solche Art von Kopplung zwischen einem fixierten (Eulerschen) 3D-Gitter und einem Lagrange'schen "Wandgitter" ursprünglich in einem frühen Versuch von Peskin 1974, die Herzdynamik numerisch zu erfassen, entstanden ist, bei dem er die sich bewegenden Wände des Herzens als interne Grenzen verwendete. Dennoch sind fast alle bisherigen CFD-Anwendungen im Bereich der Hämodynamik, mit Ausnahme der Arbeiten von Peskin und seiner Schule, weit von diesem effizienten und natürlichen Ansatz entfernt geblieben. Derartige Ansätze bringen mit sich ein hohes Potential an Vereinfachung, Genauigkeitsverbesserung und Beschleunigung der Rechenverfahren, besonders bei strukturierten Rechengittern, da die Kopplung der beiden gekoppelten Gitter einfacher zu berechnen ist.

Der oben erwähnte 3D-Euler-2D-Lagrange-Ansatz wurde unter dem Namen "immersed boundary method" technisch vielseitig entwickelt. Insbesondere wurde die ursprüngliche theoretische Einschränkung der Genauigkeit auf die erste Ordnung durch Verbesserungen am Algorithmus aus praktischer Sicht gelöst, so daß das gesamte Verfahren eine effektive Konvergenz "zweiter Ordnung" erreicht (Lai & Peskin (1999)). Ein systematischer Ansatz zur Kopplung der zwei Gitterarten (fixiertes 3D-Rechengitter für den Strömungslöser und 2D-Gitter für bewegte oder unbewegte Oberflächen komplexer Form, was auch auf Blutgefäß-Wände angewendet werden kann) wurde in den letzten Jahren unter dem Namen "immersed interface method" entwickelt und getestet. Es beruht auf einer möglichst genauen und scharfen Darstellung der Grenze zwischen dem Fluid und dem umgebenden Medium. Dagegen wird bei der Immersed-Boundary-Methode die Grenze absichtlich und kontrolliert auf eine "dünne Schicht verschmiert". Algorithmen beider Arten sind im Zusammenhang mit verschieden Diskretisationsverfahren (Finite-Differenzen, FVM, FEM) für das 3D-Rechengitter entwickelt und erprobt worden, nicht aber mit der Lattice-Boltzmann-Methode. Ihre Stabilität und Effizienz, insbesondere auch auf Parallelrechnern, ist in der Literatur bereits demonstriert worden.

Die Bestimmung der Größen, die es aus den durch die Simulation berechneten Druck- und Geschwindigkeitsfeldern zu berechnen und anschließend zu visualisieren und auszuwerten gilt, ist meistens einfach und schnell, obwohl es bei der Berechnung von räumlichen Ableitungen (die z.B. bei der Darstellung von Spannungen und Wirbel gebraucht würden) bei nicht ausreichender Auflösung zu Instabilitäten kommen kann. Um diese Probleme zu unterdrücken, werden zeitliche und räumliche Mittelung oder Glättung (Filterung) verwendet. Obwohl es eine Vielzahl von Visualisierungs-Softwarewerkzeugen gibt, die es erlauben, 3D-Strömungen darzustellen, steht nur selten, mindestens was kommerzielle Simulations-Software angeht, auch die Möglichkeit zur Verfügung, Glättungsoperationen per Mausklick durchzuführen.

Schwieriger ist die Übertragung und mehrfache Visualisierung (etwa aus verschiedenen Sichtwinkeln) der großen Datenmengen, die bei DNS entstehen. Wie bereits bei der Visualisierung von Turbulenz oft praktiziert wird, kann der Aufwand bei der Visualisierung auf ein technisch sinnvolles Niveau reduziert werden, indem die bereits berechneten Größen vor der Visualisierung reduziert (projiziert, gemittelt, strobiert, etc.) werden. In der Literatur über Hämodynamik-CFD sind derartige Überlegungen leider noch nicht bekannt.

Der Erfindung liegt die Aufgabe zugrunde, ein Verfahren zur Darstellung von dreidimensionalen, zeitabhängigen oder auch Zeit-gemittelten medizinisch verwertbaren Daten anzugeben, das detaillierte Analysen des Blutstromverlaufs auch bei Gefäßanomalien sowie unter Berücksichtigung der Pulsationen – sowohl des gesamten Blutstroms wie auch ggf. der betroffenen Gefäßteile selbst – ermöglicht und auch bei komplexer Geometrie effizient durchführbar ist. Diese Aufgabe ist durch die Erfindung bei einem Verfahren mit den Merkmalen des Anspruchs 1 gelöst. Vorteilhafte Verfahrensvarianten sind Gegenstand der Unteransprüche.

Bei dem erfindungsgemäßen Verfahren zur Darstellung von dreidimensionalen Strömungsdaten in fluidführenden Gefäßen unter Verwendung der Lattice-Boltzmann-Methode werden somit dreidimensionale geometrische, unter Verwendung eines regelmäßigen, isotropen Gitters erzeugte Strukturdaten (Geometrie-Voxeldaten) des darzustellenden Gefäßbereichs einschließlich dessen Wandung verwendet. Dieselben Datengitter oder daraus vollautmatisch erzeugte regelmäßige isotrope Rechengitter (Voxel-Gitter) werden für die Gefäßdarstellung und für die Bestimmung des Strömungsverlaufs ausgewählt werden und der Verlauf der Fluidströmung wird in zeitlichen Schritten bestimmt. Ein oder mehrere zu untersuchende Gefäßabschnitte werden ausgewählt, weiter wird ein Bereich des Gefäßes aus den dreidimensionalen geometrischen Daten ausgewählt, der den oder die Gefäßabschnitte enthält und in dem die Strömungszustände dargestellt werden sollen. Die dreidimensionalen geometrischen Daten zur Darstellung der Gefäßwände werden segmentiert, die Ein- und Ausströmrandflächen des Gefäßbereiches werden festgelegt und die Gitterrastergröße des Rechengitters wird abhängig von den Abmessungen des ausgewählten Gefäßbereiches und den Parametern der darin stattfindenden Strömung ausgewählt. Das Rechengitter wird zerlegt und die Teile des Gitters werden eliminiert, die keine durch- oder angeströmten Voxel enthalten. Die zeitliche Bestimmung der Fluidströmung wird so lange durchgeführt, bis die Strömungsdarstellung ausreichend genau ist. Die für die Darstellung vorgesehenen, d.h. interssierenden Charakteristiken der Strömung und der mechanischen Belastung der Gefäßwände werden aus der Strömungslösung bestimmt und die Ergebnisse in für die Bildgebung oder Archivierung optimierte Datenstrukturen werden übertragen und gespeichert.

Die Schnittmenge des Gefäßvolumens mit dem Rand des gewählten Bereichs wird in Ein- und Ausströmrandflächen aufgeteilt und für jeden dazu gehörenden einzelnen, ununterbrochenen Flächenschnitt werden entsprechende Arten von Randbedingungen festgelegt. Diese benötigen zeitlich und räumlich aufgelöste Randdaten. Entsprechende Angaben können aus individuellen Messungen am Patienten errechnet oder auch bereits vorhandenen Datenbanken entnommen werden. Es gibt eine Vielzahl von publizierten Verfahren zur Ermittlung, signaltechnischen und statistischen Nachbearbeitung, Modellierung und Anwendung von solchen Daten. Es werden auch geeignete Volumenkräfte angenommen, die den pulsierenden Charakter der untersuchten Blutströmungen bestimmen und auf physiologischen Messungen basiert sind.

Die Gitterrastergröße des Rechengitters wird abhängig von den Abmessungen des ausgewählten Gefäßbereiches ausgewählt und ist nicht zwangsläufig an die Größe der Geometriedaten gebunden. Bei dem erfindungsgemäßen Verfahren sind die Möglichkeit der lokalen Verfeinerung sowie auch der Vergröberung des Gitter schrittes auf der Basis einer automatischen Gebietszerlegung und Geometrieanalyse sowie einer Angabe der voraussichtlichen Flußmengen vorgesehen.

Für die zeitliche Bestimmung der Fluidströmung ist bei einigen Anwendungen eine stationäre (zeitunabhängige) Lösung der Strömungsaufgabe ausreichend, um diagnostische Aussagen zu beeinflussen. Eine solche Lösung ist einfacher und schneller als eine zeitperiodische (pulsierende) Lösung auszuführen. Die komplexe Berandung des Strömungsgebiets beeinflußt maßgeblich die gesamte Druck- und Massenflußverteilung und wird so in allen Fällen als erstes bestimmt. Zur Beschleunigung dieser Simulationsphase werden eine Linearisierung der Dynamik (Stokessche anstatt Navier-Stokes-Gleichung) und, falls sinnvoll, ein für LBM spezialisiertes Mehrgitter-Verfahren (Multigrid) verwendet. Falls erwünscht, wird daraus dann eine zeitabhängige Lösung berechnet, die durch angegebene, synchronisiert zeit-periodische Druckdifferenzen an den Ein- und Ausflußrändern und Volumenkräften im Inneren des Strömungsgebiets angetrieben wird. In dieser Simulationsphase haben die numerischen Iterationen die physikalische Bedeutung von Zeitschritten. Für LBM sind 100 bis 800 Schritte pro Pulsperiode je nach Genauigkeitsanforderungen ausreichend. Für den Kreislauf eines erwachsenen Menschen ist die Durchschnittsdauer einer Pulsperiode ca. 0,8 Sekunden, so daß jeder (explizite) Zeitschritt des LBM-Lösers etwa 1 bis 10 Millisekunden entspricht. Die Konvergenz aller interessierenden Größen auf einen zeitperiodischen Verlauf wird stets gemessen und die Simulation dann unterbrochen, wenn eine ausreichende Näherung an ein periodisches Verhalten erhalten wird oder ein Zeitlimit überschritten wird. Noch bei Simulationsanfang wird dem Benutzer automatisch angekündigt, welche Anzahl von simulierten Pulsperioden diesem Zeitlimit entspricht und ob diese nach vorhandenen Regeln und Heuristiken ausreicht, um die gewünschten Größen mit der gewünschten Genauigkeit berechnet zu können. So wird eine Optimierung des Rechenaufwands bereits vor der möglicherweise aufwendigen Simulation ermöglicht.

Das erfindungsgemäße Verfahren kann gemäß verschiedenen Alternativen entweder schnell oder hochgenau durchgeführt werden. Es können beispielsweise bewegte Wände simuliert werden, wobei das Verfahren genauer, aber wesentlich aufwendiger als bei der Annahme von starren Wänden ist. Bei kleineren und besonders bei intracraniellen (sich innerhalb der Schädelhöhle befindenden) Gefäßen sowie beispielsweise, wenn die Beschaffenheit der Gefäßwände unbekannt ist, kann auf die Mitberechnung der Wandbewegung verzichtet werden.

Ein alternatives Einsatzszenario sieht vor, die Gefäßbewegung, die einschließlich von, aber nicht nur durch die Wandpulsation bestimmt ist, aus mehreren und somit Zeit-aufgelösten Tomographieaufnahmen zu ermitteln, anstatt sie als Teil der modellierten und simulierten Dynamik bei jedem Zeitschritt der Simulation neu berechnen zu müssen. Bei diesem Szenario ist also die Kopplung zwischen Fluid und Wand einseitig: Die Wandbewegung induziert eine Strömung, aus der die Wandbelastung dann ausgerechnet werden kann, diese Belastung wird aber dann nicht als Antrieb der weiteren Wandbewegung betrachtet und es wird keine Wandbewegung simuliert, da sie aus den Aufnahmesequenzen bereits bekannt ist.

Wenn Wand-Scherspannungen bestimmt werden sollen, müssen die Wände mit hoher Genaugikeit dargestellt werden, wofür LBM-Randbedingungen zweiter Ordnung, die aus der Literatur bekannt sind, verwendet werden können. Diese und weitere, alternative Rechenverfahren im Wandbereich sind aufwendiger als das einfache Bounce-back-Verfahren (BBL) und führen im Gegensatz zu ihm zu einem geringen "Massenverlust" wärend der Simulation. Doch bei ausreichender Auflösung sind weder dieser Verlust noch der Fehler bei der diskretisierten Variante der Inkompressibilität der Strömung bedenklich. Bei gewissen Fragestellungen, z.B. Behandlungen von intracraniellen Aneurysmen oder bei Analysen der Durchströmbarkeit von ganzen Zweigen des Gefäßsystems, kann auf das BBL zurückgegriffen werden.

Meistens werden die dargestellten Gefäßabschnitte einen Gefäßstrukturteil enthalten. Hierunter sollen insbesondere Deformationen und Störstellen verstanden werden, beispielsweise Implantate oder Geometrieänderungen wie Aneurysmen und Stenosen. Diese können räumlich vorhanden sein oder anhand von gespeicherten Datensätzen synthesiert und in die Gefäßdarstellung eingebaut werden, um die sich entsprechend ergebenden Strömungs- und/oder Druckverhältnisse zu bestimmen. Umgekehrt können auch in der Realität vorhandene Gefäßdeformationen virtuell von der Gefäßdarstellung entfernt werden, um beispielsweise die Effizienz und die langfristigen Auswirkungen einer Erweiterung des durchströmbaren Lumens zu untersuchen. Letzteres kommt in Frage z.B. bei Kranzgefäßoperationen zur Linderung von Stenosen oder bei "Clipping" von Aneurysmen und von arteriellen Zweigen in der Neurochirurgie.

Die Verwendung ausschließlich von regelmäßigen 3D-Datengittern bei dem erfindungsgemäßen Verfahren führt zu einem Sprung in der Effizienz und der Standardisierung von simulationsgestützten Untersuchungen. Dies basiert auf den folgenden Merkmalen: Es werden dieselben und als ausreichend feinmaschig angenommenen Rastergitter direkt oder nach einer vollautomatischen und sehr schnellen Abbildung bzw. Interpolation in der Koordinatenrichtung, die der Hauptache des Aufnahmegeräts entspricht, auf kollineare, teilweise überlappende Gitter mit kubischen Voxeln verwendet. Diese Voxelgitter können dann für die Bestimmung des Strömungsverlaufs unmittelbar verwendet werden. Eine Spezifikation der Abbildung in Abhängigkeit von der Größe, Krümmung und Verzweigung der betrachteten Gefäßstruktur ist Bestandteil des Verfahrens. Es ist keine aufwendige geometrieangepaßte Generierung von unstrukturierten 3D-Rechengittern erforderlich, um alle Einzelheiten der Gefäßgeometrie mit guter Auflösung darzustellen. Dies hat dramatische Verminderungen der Komplexität und der Gesamtdauer der Simulation zur Folge.

Die Rekonstruktion der Gefäßgeometrie bezieht sich auf mehrere Raumteile, die sich in ihrem Abstand von der untersuchten Gefäßmalformation sowie auch in der jeweils benötigten Auflösung und Modellierungspräzision unterscheiden. Die Rekonstruktion ist schnell auf aktuellen Rechnern, kann aber nur teilweise automatisiert werden. Die Abmessung des zu untersuchenden Gefäßabschnitts, z.B. der Umfang einer Malformation selbst und der Raumteile von klinischem, Interesse in dessen Umgebung ist von individuellen Unterschieden in der Gefäßgeometrie und der Art und Ursachenlage der Beschwerden mitbestimmt und kann interaktiv bestimmt werden.

Um eine robuste Strömungssimulation mit der beschriebenen Art von regelmäßigen Gittern bei allen zu erwartenden Geometrien zu ermöglichen, wird eine angepaßte Variante der Lattice-Boltzmann-Methode zum Einsatz gebracht, bei der das Rechengitter zerlegt wird und die Teile des Gitters eliminiert werden, die keine durch- oder angeströmten Voxel enthalten. Im Vergleich mit CFD-Standardverfahren wird dadurch einerseits eine hohe Effizienz auch auf hochparallelen Rechnerarchitekturen und andererseits weniger Abhängigkeit von der Komplexität der Geometrie und des Rechengitters erzielt.

Der gesamte Arbeitsvorgang kann anhand von Meßdaten vollzogen werden, die ausschließlich mit minimal-invasiven und nicht-invasiven Verfahren erfaßbar sind.

Das Verfahren liefert vierdimensional (zeitlich und 3D räumlich) aufgelöste Vorhersagen über den Verlauf der Blutströmung durch Arterien in der Nähe von typischen Malformationen wie Stenosen oder Aneurysmen. Dadurch können auftretende mechanische Belastungen und der Zusammenhang zwischen Gefäßgeometrie, Hämodynamik und Wandbelastung, die durch physiologische Wandmodellierung wiederum die Gefäßgeometrie beeinflußt, untersucht und besser vorhergesagt werden.

Das erfindungsgemäße Verfahren kann bei der individualisierten Untersuchung möglicher mechanischer Ursachen, die mit der ge nauen Raumstruktur der Blutgefäße eines Patienten verbunden sind, für die Entstehung von Gefäßmalformationen und dessen Weiterentwicklung angewendet werden. Aufgrund der hohen Effizienz können so Diagnose und Prognose individuell für mehrere Patienten an einem Tag mit hämodynamischen Daten unterstützt werden, die wesentlich umfangreicher und detaillierter im Vergleich zu direkten Messungen sind und darüber hinaus keine zusätzliche Manipulationen am Patienten erfordern. Es ist ein regulärer klinischer Einsatz der Strömungsuntersuchungen und auch nachträgliche Kontrollen von Untersuchungsdaten möglich.

Die direkte Verwendung von Tomographiedaten im regelmäßigen Voxelformat erlaubt eine einfache und gleichzeitig zuverlässige Untersuchung der Wirkung bereits im Körper befindlicher Implantate (Stents, Coils, usw). Da Modelle für die biomechanische Reaktion von Blut (etwa durch Koagulation oder Verdünnung) und von Gefäßwänden auf mechanische Spannungen zur Verfügung stehen, kann die Wahrscheinlichkeit von klinisch wichtigen Folgen der Einführung von Implantaten durch eine Serie von Simulationen individuell eingeschätzt werden. Beispiele solcher Folgen sind die Hyperplasie und Restenose, die Wochen nach einer Stentimplantation auftreten können, die erwünschte Verdichtung eines Aneurysmenlumens infolge einer durch Implantate (z.B. GDC) verursachten Blutkoagulation oder die später mögliche unerwünschte Rekanalisierung des Aneurysmenlumens, sowie in jedem Fall eine Gefahr von Embolien stromabwärts von dem Ort des Implantats.

Eine weitere wichtige Anwendung ist die "virtuelle Angioplastie", die eine detaillierte Vorhersage der klinischen Folgen von geplanten, aber noch nicht operativ durchgeführten Einsätzen von Implantaten oder von vasoplastischer Chirurgie erlauben wird. Die tomographisch aufgenommene Gefäßgeometrie kann nach softwaregestützter Erkennung der räumlichen Struktur der Gefäße virtuell verändert werden, z.B. indem Gefäße verschlossen oder mit Implantaten versehen werden. Die entstehende neue Gefäßgeometrie entspricht dann nicht dem aktuellen Zustand des Patienten, sondern der Vorstellung des Arztes, wie die Blutgefäße des Patienten nach einem Eingriff aussehen sollen. Diese virtuelle Gefäßgeometrie kann als Grundlage für numerische Simulationen und auf dessen Basis für quantitative Aussagen über die Hämodynamik und Wandbelastung dienen, genau so gut, wie die ursprüngliche, direkt aus Tomographiedaten gewonnene Geometrieinformation.

Die Einschränkungen im Umfang möglicher Anwendungen sind im wesentlichen auf Einschränkungen bei den Aufnahme- und Registrationsverfahren zurückzuführen. Der Herz- und Lungen-Bereich unterliegt einer nicht unterdrückbaren starken Bewegung. Durch die Erfindung ist es möglich, bei Einsatz der Verfahrensmöglichkeiten, insbesondere der dynamischen Kopplung von Blutstrom und Wandbewegung, auch die Blutströmung im Herzen und/oder den Kranzgefäßen zu berechnen. Venen sind schwer zu modellieren und oft auch schwer mit Tomographie genau zu lokalisieren. Kleine Gefäße (Stenosen unter 2 mm Durchmesser und Aneurysmen unter 4 mm) sind aufgrund begrenzter CT-Auflösung (ca. 0,5 mm) zur Zeit nicht zuverlässig abschätzbar.

Die Strömungssimulation ist derzeit das zeitaufwendigste Element des Verfahrens. Mit dem Einsatz der Lattice-Boltzmann-Methode und von kleineren Hochleistungsrechnern (2 bis 8 Prozessoren), die nur für einen Bruchteil der Kosten eines Tomographen (z.B. CT oder MRI) zu erwerben sind, kann die Dauer einer typischen Simulation in einfacheren neurochirurgischen Anwendungen bereits auf 1 Stunde reduziert werden. Die Vorbereitung der Simulation einschließlich Segmentierung, 3D-Rechengitter, Verfeinerung, Angabe von Randbedingungen, usw. soll im Durchschnitt 5 bis 10 Minuten betragen und die Auswertung der Simulationsergebnisse 5 bis 20 Minuten, incl. Archivierung der benutzten Gefäßgeometrie, Randdaten, des konvergierten Strömungszustandes, der für weitere Simulationen, z.B. für eine "virtuelle Angioplastik" gebraucht werden kann, sowie ausgewählte, über die Zeit integrierte Simulationsergebnisse. Das Verfahren ist damit nicht für Notfallbehandlungen geeignet, dagegen aber bei der Therapieplanung und der nachträglichen Kontrolle als Standardmittel zur klinischen Unterstützung verwendbar.

Die Methode kann auch als Upgrade für bereits im Einsatz befindliche Tomographieapparaturen eingebaut werden. Benötigt werden dazu aktuelle leistungsstarke Rechnercluster und herstellereigene Interfacesoftware, die eine Verbindung mit der Meßapparatur des Herstellers ermöglicht. Für spezifische Anwendungsbereiche wie z.B. Neurochirurgie ist für die zuverlässige Simulation eine Scannerauflösung im Submillimeterbereich notwendig, die bei manchen Tomographiegeräten nicht erreichbar ist.

Der Hauptvorteil des erfindungsgemäßen Verfahrens für eine simulationsgestützten Therapieplanung und Kontrolle besteht in der Automatisierung, Standardisierung und Beschleunigung von Blutstromsimulationen in Adern. Eine Segmentierung der dreidimensionalen geometrischen Daten ermöglicht die Darstellung der Gefäßwände. Von entscheidender Bedeutung ist der Einsatz von regelmäßigen Datengittern bei allen Stufen der rechnerischen Bearbeitung und insbesondere bei der Lattice-Boltzmann-Methode (LBM) für die Strömungssimulation. Alle bisher publizierten Simulationen, die auf individuellen Tomographiedaten von Patienten basieren, sind mit anderen Verfahren (FVM, FEM) durchgeführt worden. Es sind nur einzelne Versuche bekannt, LBM zur Blutstromsimulation einzusetzen, der Aufbau der durchströmten Geometrie ist dabei stets nicht in Zusammenhang mit Daten von individuellen Patienten gebracht worden. Auch bei den Simulationen mit herkömmlichen Verfahren (FVM, FEM) ist der Einfluß des Strömungsgebietsumfangs nicht untersucht worden. Der Einfluß von Randbedingungen am Ein- und Ausströmort und von der räumlichen Auflösung ist nur selten und dann nur empirisch untersucht worden, ohne daß dabei eine Zuverlässigkeitsregel (z.B. in Form eines Fehlerintervalls) festgelegt worden wäre. Die virtuelle Chirurgie ist als zukunftsträchtiges Mittel der simulationsgestützten Therapieplanung bereits diskutiert worden, allerdings nur in einem Zusammenhang mit Koronararterien und Aortensegmenten. Dabei ist die gegenwärtige Machbarkeit von, 3D-Simulationen für die klinische Praxis negativ eingeschätzt worden.

Erst das erfindungsgemäße Verfahren mit räumlich dreidimensional aufgelöster Blutstromsimulation stellt eine mit vorhandener Hardware umsetzbare Option für simulationsgestützte Diagnostik und Therapieplanung von Gefäßmalformationen dar. Wichtig ist die Festlegung des erforderlichen Umfangs (eine ausreichend weite Umgebung der Malformation, je nach Typ, nach Komplexität der Gefäßstruktur in der Umgebung, und nach individuellen Besonderheiten) und die erforderliche (räumliche und zeitliche) Auflösung der Blutstromsimulation und dementsprechend die der mit Tomographieverfahren gewonnenen Ausgangsdaten. Des weiteren entscheidend ist die durchgehende Verwendung nur von regelmäßigen Datengittern, was sowohl die Effizienz als auch die Genauigkeit bei den einzelnen Schritten des Verfahrens verbessert.

Für die Durchführung des Verfahrens werden keine Fachleute auf dem Gebiet der Strömungsmechanik mehr gebraucht. Es werden Schnittstellen für die Einführung von biomechanischen, zytologischen und anderen (für die Darstellung der normalen und krankhaften Remodellierung von Gefäßen relevanten) Modellen vorgesehen. Dadurch können zusätzliche wichtige Einflüsse auf die Hämodynamik und Gefäßwandbelastung, z.B. durch Wandpulsation, Blutrheologie und Koagulation, transiente und systematische Änderungen in dem biochemischen Verhalten von Wandzellen usw. mitberechnet werden.

Das erfindungsgemäße Verfahren eignet sich als Standardmethode zur simulationsgestützten Diagnostik, Therapie und Kontrolle von Blutgefäßmalformationen, somit zur Diagnostik, Therapieplanung und Kontrolle von vielen Symptomen in der Kardiologie, Neurologie und Gefäßchirurgie algemein und ist im Gegensatz zu den bisherigen 3D-Blutstromsimulationsverfahren in jeder aktuell ausgerüsteten Klinik verwendbar.

Das erfindungsgemäße Verfahren basiert auf folgenden Vorausset zungen:
Es müssen räumlich aufgelöste Tomographiedaten (3D von näherungsweise unbewegten Körperteilen oder 4D, also einschließlich Zeitabhängigkeit, bei Herz und Lungen) vorhanden sein. Die Möglichkeiten der vorhandenen Tomographiehardware beschränken dadurch die Anzahl und Lage der behandelbaren Gefäße.

An allen proximalen Endquerschnitten (in bezug auf die vorhandene räumlich beschränkte Abbildung, die als Geometrieangabe für die numerische Simulation benutzt werden soll) von Arterien muß eine unabhängige zeitlich aufgelöste Meßinformation über den Volumenstrom oder den Blutdruck vorliegen. Blutflußinformation kann z.B. durch Ultraschallmessungen mit Doppleranalyse (ein schnelles Verfahren, das bei ausreichend großen Gefäßen standardmäßig eingesetzt wird) oder durch MRI-Doppler-Verfahren gewonnen werden. Zur Abmessung des Druckverlaufs in größeren Adern steht eine Reihe von meistens nichtinvasiven Verfahren zur Verfügung.

Modelle für das mechanische Verhalten von gesunden und beschädigten Gefäßwänden, sowie für den räumlichen Verlauf der Lumen-Querschnittsfläche kleinerer Gefäße werden für eine vollständige Modellierung und Simulation benötigt. Entsprechende Information ist in der Fachliteratur vorhanden, siehe z.B. G. A. Holzapfel and R. W. Ogden, eds., "Biomechanics of soft tissue in cardiovascular systems", CISM courses and lectures 441, Springer, Wien, 2003. Für manche Anwendungen ist sie allerdings noch nicht ausreichend und so werden dann entweder 4D-Aufnahmen oder die Annahme von starren Wänden in Zusammenhang mit 3D-Aufnahmen verwendet.

Im folgenden wird das erfindunsgemäße Verfahren mehr im einzelnen erläutert. Der Kern des Verfahrens ist eine direkte numerische Simulation (DNS) des Blutstroms in der Nähe von arteriellen Malformationen, die mit Standardverfahren der Angiographie individuell und zuverlässig bei Patienten festgestellt werden können. Das Verfahren beruht auf drei Hauptschritten, die im wesentlichen bei jeder Anwendung von CFD-Software durchgeführt werden müssen.

Pre-Processing: Die Geometrie des durchströmten Raumgebietes (Gefäßsystem) und die passenden Randbedingungen werden festgelegt. Entsprechend wird ein System von Rechengittern ausgelegt, die zusammen das gesamte durchströmte Volumen decken und die erforderliche räumliche Auflösung lokal gewährleisten, wobei jedes einzelne Gitter nur kubische Voxel einer entsprechenden Größe hat und die minimale und maximale Gittergröße an die Architektur des verwendeten Rechners bei der Installation der Software optimal angepaßt werden kann.

Processing: Eine numerische Simulation der pulsierenden Blutströmung durch die festgelegte Geometrie wird so lange durchgeführt, bis eine verläßliche Auswertung der Strömungsparameter von Interesse möglich ist (numerisch konvergierte Lösung im weiter oben beschriebenen Sinne vorhanden ist), einschließlich raum- oder möglicherweise auch zeitgemittelte Statistiken.

Post-Processing: Die Simulationsdaten werden ausgewertet und visualisiert. Räumlich verteilte Felder, z.B. Schwankungsamplituden der Scherspannungen an gewählten Abschnitten der Gefäßwände oder "Turbulenzintensität" nach Zeit-Mittelung im Inneren eines Gefäßlumens, sowohl wie auch einfache Zeitreihen, z.B. Des Druckverlaufs an Kontrollpunkten innerhalb von Aneurysmen, und charakteristische Einzelwerte, wie z.B. räumlich gemittelte Oszillationsindizes oder Massenfluß-Verteilungen an Gefäßbifurkationen, werden berechnet und für die Abschätzung des Einflusses von hämodynamischen Faktoren, aber auch von Ungenauigkeiten bei der Modellierung und der numerischen Approximation der berechneten Dynamik bereitgestellt.

Bei jedem dieser drei Hauptschritte im Arbeitsprozeß laufen folgende neue Maßnahmen ab, die die Effizienz und Genauigkeit erhöhen und die Standardisierung erleichtern: Die Verwaltung aller Datengitter ist vollautomatisiert. Alle dreidimensionalen (3D) Daten werden ausschließlich auf regelmäßigen, rechteckigen und so oft wie möglich auf isotropen Gittern dargestellt. Dadurch können die schnellsten und höchstpräzisen Versionen der jeweiligen Algorithmen zur numerischen Glättung, Differenzierung, Integrierung, usw. eingesetzt werden.

Für die DNS der Navier-Stokesschen Strömungsmechanik wird eine Lattice-Boltzmann-Methode eingesetzt. Dadurch wird auch bei sehr komplizierten Gefäßgeometrien einerseits eine optimale Rechenleistung auch auf Parallelrechnern gesichert und andererseits eine zuverlässige Bestimmung aller hämodynamischen Parameter gewährleistet, insbesondere des Drucks in vielfach verzweigten Gefäßsystemen und der viskosen Spannungen in der Nähe von Gefäßwänden. So ist aus eigener Forschung bekannt, daß es zur DNS Auflösung auch der am stärksten fluktuierenden Moden (Womersleyzahlen bis ca. 20), die zur realistischen Darstellung einer pulsierenden Strömung besonders in Wandnähe noch wichtig sind, eine mäßige räumliche Auflösung von 15 bis 25 Gitterpunkten im Gefäßdurchmesser bei Verwendung von LBM ausreicht, wobei die Wandschichten, in denen sich die höchsten Moden lokalisieren, mit nur 3 bis 5 Gitterpunkten dargestellt werden. Eine entsprechende Gittergröße in der "Region of Interest" kann (automatisch) und soll während des Preprocessing gewährleistet werden und kann danach während der Lösungsschritts leicht von der LBM bearbeitet werden.

Es sind Vergleiche zwischen LBM und FVM bekannt, die wesentliche Vorteile der LBM erkennen lassen, z.B. Diese zeigen, daß bei Strömungsgebieten mit komplexer Geometrie die LBM eine wesentlich geringere Anzahl von Gitterpunkten braucht, um eine vergleichbare Genauigkeit bei der Ermittlung des stationären Druckverlusts zu erreichen. Bei turbulenten Strömungen hat ein Vergleich mit Pseudospektralverfahren, die von höchster Genauigkeit sind, eine vergleichbare Genauigkeit der berechneten Statistiken sowie Geschwindigkeitsvorteile bei Gittergrößen über etwa 100 Punkte in einer Koordinatenrichtung. Da etliche Spektralverfahren nur sehr schwer bei komplexen Geometrien einsetzbar sind, konkurriert LBM in solchen Fällen mit FEM und FVM. Aus der Literatur ist es bekannt, daß FEM sehr zeitaufwendig und recht kompliziert in der Anwendung auf Blutgefäßströmungen sind. Eigene Vergleiche mit kommerzieller CFD-Software auf FVM Basis haben gezeigt, daß die Berechnung einer stationären Lösung mit der FVM-Software, die Gebrauch von fortgeschrittenen Mehrgitter-Verfahren (Multigrid, MG) macht, eine mit LBM ohne MG vergleichbare Rechenzeit in Anspruch nimmt. Durch den Einsatz von MG für LBM kann bis zu eine Größenordnung an Beschleunigung erreicht werden. Die grundsätzliche Beschreibung des LBM-MG ist bereits im Internet zu finden zusammen mit einer begründeten Auswahl einer bestimmten Vorgehensweise aus einer Reihe von sich theoretisch anbietenden Alternativen. Bei der Simulation von pulsierenden Strömungen hat derselbe Vergleich einen Vorsprung der LBM (ohne MG) von wenigstens einer Größenordnung in der Rechenzeit bereits gezeigt. Mittels eines solchen Mehrgitter-Verfahrens können die stationären Skalar- und Vektorgrößen für die Lattice-Boltzmann-Methode schneller bestimmt werden.

Durch diese und weitere unten beschriebene Merkmale ist das erfindungsgemäße Verfahren viel einfacher, weitgehend automatisierbar und wesentlich präziser und schneller bei der Ermittlung von klinisch relevanten Strömungsmustern und hämodynamischen Belastungen.

Für das Pre-Processing werden zwei Arten von Eingaben von dem hier beschriebenen Verfahren benötigt: 3D-Tomographie-Daten, die das Blutgefäßgewebe unterscheiden lassen, und zeitabhängige Messungen des Blutstroms in Adern, die sich proximal in der Nähe der untersuchten Malformation befinden. Um eine zuverlässige Geometrieerkennung und danach auch Strömungssimulation zu gewährleisten, werden Anforderungen an die Auflösung der Eingangsdaten gestellt: Ein Pulszyklus muß durch die Angabe von mindestens 16 unabhängigen Datenaufnahmen, die unterschiedli chen Zeitphasen in einem Zyklus entsprechen, davon mindestens 8 in der Systole und 8 in der Diastole bei einem gesunden Herzzyklus, mit zureichender Zeit-Auflösung definiert sein. Die räumliche Auflösung von Geschwindigkeitsprofilen darf nicht mehr als 3 mal gröber (das Verhältnis bei einem Vergröberungs/Verfeinerungs-übergang bei benachbarten Gitterblöcken) sein als die des Rechengitters, das die Rand-Daten empfangen soll. Im Gegenteil muß auf die bei CFD-Software übliche globale Randbedingungsart, bei der ein Gesamtmassenfluß durch einen gewählten Querschnitt angegeben wird, zurückgegriffen werden.

Die zur Erstellung von Eingabedaten typischerweise verwendeten Tomographie-Daten (CT, MRT, etc.) werden als 3D-Datensatz aus Grauwerten, die auf einem regelmäßigen Gitter im physikalischen Raum (und nicht in dem Frequenzraum) gegeben sind, erwartet, in Zusammenhang mit entsprechenden Grauwertintervallen und Filterfunktionen, die wie bei der in der medizinischen Bildgebung üblichen Segmentierung zur Angabe der Unterscheidungsschwellen benutzt werden. Die Vorbereitungsphase, in der diese Angaben vom Anwender erst festgelegt werden müssen, kann so mit bestehender, herstellerspezifischer medizinischer Software durchgeführt werden. Üblicherweise liegen in der Klinik solche Grauwertdaten in einer Variante des DICOM-Formats vor, typischerweise als eine Vielzahl von 2D-Aufnahmen in parallelen Schichten vorgegeben. Die nominale Auflösung (der aus der DICOM-Datei abzulesende räumliche Gitterschritt) innerhalb einer Schicht darf dabei in einem Verhältnis der Gitterschrittgrößen von nicht weniger als 1:3 zur Auflösung in der dritten, axialen Koordinatenrichtung stehen. (Das ideale und bisher maximale Verhältniss von 1:1 wird erst in den neuesten Geräten und dann nicht bei allen Meßverfahren erreicht.) Die Überprüfung, ob die so angegebene Auflösung der tatsächlichen technischen Auflösung während der Datenaufnahme entspricht, kann nicht immer automatisch erfolgen und muß gegebenenfalls vom Anwender durchgeführt werden. Es wird also angenommen, daß die axiale Auflösung nicht wesentlich schlechter ist als die entlang den Schnittebenen. Die gröbere von den beiden Auflösungen muß aus reichend sein, um alle Durchmesser von Gefäßen, die in dem Gebiet von Interesse (Region of Interest) mitberechnet werden sollen (also in der segmentierten Geometriebeschreibung) mit wenigstens 4 Voxeln in der Eingabedatei erscheinen zu lassen. Man hat z.B. eine Ader von 10 mm Durchmesser, die einen Abschnitt mit 75 % Stenose vorweist. Wenn dieser Abschnitt im interessierenden Gebiet liegt, darf die Auflösung nicht gröber als 0,6 mm pro Pixel sein. Bei Aneurysmen muß ihr Hals mit wenigstens 4 und der kürzeste Durchmesser ihres Sacks mit wenigstens 6 Voxeln dargestellt sein. Bei einer isotropen Auflösung von 0,4 mm, wie sie zur Zeit mit dem neuen 64-Zeiler-CT-Gerät von Siemens erreichbar sein soll, können Aneurysmen bis 3 mm und Adern von 1,5 bis 2 mm behandelt werden.

Die zweite Art von erwarteten Angaben sind Puls-, Blutdruck- und Blutstromvolumen-Daten, die es erlauben, Einström-Randbedingungen für die DNS mit ausreichender Präzision festzulegen. Eine Reihe von nichtinvasiven und minimal-invasiven, weitgehend standardisierten Verfahren stehen für die unabhängige Bemessung dieser Parameter zur Verfügung. In G. Pedrizzetti and K. Perktold, eds., Cardiovascular fluid mechanics, CISM courses and 1ectures 446, Springer, Wien, 2003 kann z.B. eine Übersicht bestehender Ultraschall-Meßverfahren zur Bestimmung von Geschwindigkeiten und Scherspannungen gefunden werden. Die räumliche Verteilung der Blutstrom-Geschwindigkeit in Querschnitten von mittelgroßen und kleineren Adern, die für die Angabe von Ein- und Ausström-Randbedingungen benötigt wird, ist für standardisierte klinische Anwendungen auch in vorhersehbarer Zukunft nicht als meßbar zu betrachten, mit der möglichen Ausnahme von großen Adern, die nahe unter der Haut und nicht von hartem Gewebe gestützt liegen, was z.B. bei der Carotisarterie bis über ihre Bifurkation hinaus auch gegeben ist. Für die haufigsten Fälle von tiefer liegenden Adern muß eine Datenbank erstellt werden, die aus der lokalen Reynoldszahl (aus Ultraschall oder sine-contrast MRI oder auch Schätzwerten) und Gefäßkrümmung (aus CT, MRT) ein plausibles Geschwindigkeitsprofil errechnen läßt. Dazu müssen eine Reihe systematischer hydrodynamischer Berechnungen und anschließender medizinischer Verifikationstests durchgeführt werden. Diese Daten sollen eine räumliche Auflösung von mindestens 20 Voxeln im Lumendurchmesser vorweisen. Die optimale zeitliche Auflösung, die für die Angabe der systolischen Pulsphase in der DNS benötigt wird, entspricht einem gleichmäßigen Zeitschritt von 4% der Gesamtdauer einer Pulsperiode oder auch kürzer, bis 1%, in Bereichen mit sehr steilem oder kompliziertem Pulsverlauf, wie z.B. in den ersten Aortenverzweigungen. In den meisten praktischen Fällen ist diese Genauigkeit nicht erreichbar, besonders bei tief liegenden Adern. Auch dieses Problem soll durch eine Datenbank von Pulsmustern gelöst werden, die es erlaubt, aus einfachen Druckverlauf-Messungen (siehe weiter oben Angaben minimaler Zeitauflösungen) mit Einbezug von weiteren Angaben, wie Alter, Geschlecht und algemeinem Gesundheitszustand des Patienten, Lage des Gefäßes etc., einen plausiblen Pulsverlauf in der Nähe der untersuchten Malformation zu errechnen. Es bleibt nicht festgelegt, ob für diesen Zweck auch Korrelationen und TFM (transfer function method) für bestimmte Körperteile zuverlässig zum Einsatz gebracht werden können.

Geometriebestimmung: Um den Einfluß von unvermeidlichen Fehlern bei der Angabe von Ein- und Ausström-Randbedingungen (EARB) zu vermindern, sollen diese Randbedingungen nur weit genug von der untersuchten Malformation angegeben werden. Dabei wird damit gerechnet, daß wirbelerzeugende Merkmale der Gefäßgeometrie (wie Bifurkationen und starke Krümmungen) die lokale Struktur der Strömung stark beeinflussen und damit den Einfluß von Details bei der Angabe der EARB reduzieren. Sowohl aufwärts als auch abwärts der Hauptrichtung der Blutströmung sollten die Ein- bzw. Ausström-Randflächen weiter als die ersten zwei Krümmungen und der ersten Bifurkation von dem Ort der Malformation entfernt liegen. In den Fällen, wenn die Anatomie oder der Umfang der vorliegenden Tomographiedaten dies ausschließt (typischerweise im Einströmungbereich), wird eine Untersuchung des Einflusses der EARB auf die Strömungsstruktur der entsprechenden numerischen Lösung empfohlen, die im Vergleich zu wenig stens zwei Varianten der Simulation, die sich nur durch das Profil oder die Art (Druck oder Massenfluß) der Einströmungs-Randbedingung unterscheiden, besteht. Die automatische Verifikation der EARB Bedingungen ist möglich, erfordert allerdings eine (wiederum automatische) topologische Gefäßstruktur-Erkennung und dafür auch eine gute räumliche Auflösung der Tomographie-Daten weit distal von der Malformation. Falls eine solche Auflösung nicht erreichbar ist, oder die Gefäßstruktur außerordentlich kompliziert ist, muß die Auswahl der Randflächen von einer Bedienungsperson, z.B. einem Arzt, getroffen werden. Diese Notwendigkeit wird automatisch signalisiert. Es wird zusätzlich verlangt, daß alle größeren Gefäße (im Gegensatz zu "kleineren Gefäßen", siehe unten) mit dargestellt werden, die mit dem Ort der Deformation (Stenose, Aneurysma, Stent, Bypass, etc.) innerhalb einer Kugel mit Durchmesser 3 L direkt mit der Deformation verbunden sind, wobei L die Länge (den größten Durchmesser) der untersuchten Deformation bezeichnet. Diese Kugel muß zum großen Teil in dem interessierenden und für die Strömungsbestimmung berücksichtigten Gebiet enthalten sein, Gefäße mit einem Durchmesser unter 30% des effektiven Durchmessers der größten Ader, die an der Deformation direkt vorbeiläuft, werden als "kleinere Gefäße" bezeichnet. Solche Gefäße werden aus der Geometriebeschreibung automatisch entfernt, wenn sie nicht direkt in oder vorbei an der Deformation führen. von den verbliebenen Gefäßen "endet" die Mehrzahl an Randflächen, and denen Ein- oder Ausström-Randbedingungen vorgeschrieben werden müssen, obwohl sie unbekannt sind. An diesen Flächen werden synthetische Randbedingungen vorgegeben, die dem Massenfluß proportional zur Querschnittfläche des jeweiligen Gefäßes (in bezug auf seine Achse, die von dem oben erwähnten Erkennungsverfahren für die topologische Struktur der vorhandenen binär segmentierten Gefäßgeometrie bereits automatisch berechnet oder mit Entscheidung eines überwachenden Arztes festgelegt worden ist) und dem dynamischen Massenerhaltungsgesetzes entsprechend bestimmt werden. Wichtige Gefäße können weit Stromabwärts (distal) nach der Deformation aus physiologischen oder auch aus technischen Gründen schlecht erkennbar werden. In diesen Fällen wird eine Bedienungsperson, z.B. der Arzt, nur bei der Festlegung des Verlaufs der Gefäßachse gefragt und unterstützt, während der Gefäßdurchmesser anhand von bekannten Näherungsformeln exponentieller Art extrapoliert wird. Um einen plausiblen Verlauf der Achse solcher Gefäße zu errechnen, werden lokale Schwellwertanpassungen, Regularisierung z.B. mit nichtlinearer Diffusion, angepaßte Einbeziehung des Grauwertgradienten in die Ortung der Gefäßwand, sowie weitere bekannte fortgeschrittene Segmentierungsverfahren einbezogen.

Simulationen von Gefäßen mit Implantaten (Stents, Coils, aber auch Bypässe) werden mit dem neuen Verfahren auch individuell patientenabhängig möglich. Die Geometrie und die mechanischen Eigenschaften des Implantats sollen dafür bekannt sein und eine ausreichend aufgelöste 3D-Tomographie-Aufnahme mit dem bereits plazierten, deutlich erkennbaren Implantat vorhanden sein. Wie auch bei der Darstellung der Wände (als gekrümmte Oberflächen ohne Dicke) und der Achsen (als gekrümmte Linien im 3D-Raum) von modellierten Blutgefäßen, wird die eigentliche Lage des Implantats durch gekrümmte geometrische Objekte ohne Dicke "registriert" (im Jargon der medizinischen Bildgebung). Diese werden als Deformationen entsprechender, voraus vorbereiteter geometrischer Modelle des Implantats in einer Referenzlage (z.B. ein Stent, der in einem geraden oder einfach (mit konstantem Krümmungsradius) gebogenen Rohr mit konstantem, kreisförmigem Querschnitt disponiert ist), die dem vorliegenden Tomographie-Datensatz in einer geeigneten Norm am nächsten kommen, erst automatisch bestimmt und danach mit Unterstützung von bildgebenden Verfahren und durch Simulation ihrer mechanischem (hyperelastischen, plastischen, oder auch, z.B. bei natürlichen Implantaten, viskoelastischen) Bewegungen zwischen den Gefäßwänden, dessen Rheologie wie bei der Berechnung der dynamischen Kopplung von Wand und Strömung modelliert wird. Nachdem die aktuelle räumliche Lage des Implantats bestimmt ist, können anhand seiner bekannten mechanischen Eigenschaften oder auch durch einfache "Verkrümmungsexperimente" mit räumlichen Abwei chungen, die von der Erkennungssoftware angegeben werden, die mechanischen Spannungen innerhalb des Implantats sowie auch zwischen ihm und den Gefäßwänden ermittelt werden.

Der fehlende, hämodynamische Anteil dieser Belastungen wird dann mittels DNS eingeschätzt. Dabei werden, bei ausreichender Auflösung der DNS, Coils mit einfacher Geometrie im aktuellen Zustand sowie Stents als 3D-Körper dargestellt, in dem ihre angepaßten "Skelette" auf eine Breite ausgebreitet ist, die der eigentlichen Dicke der Drahtquerschnitte möglichst genau entspricht, aber nicht dünner als 3 Voxel auf dem Rechengitter ist. Bei geringer Auflösung der DNS müssen Coil-Bündel als poröses Medium modelliert werden, wobei die Porosität aus den Tomographie-Aufnahmen (bei bereits im Körper vorliegenden Implantaten) oder aus eigener Erfahrung einer Bedienungsperson, z.B. des Arztes (bei der medizinischen Planung), geschätzt werden müssen.

Feinere Coils sowie Gele, Verklottungen, weiche Plaquen und ähnliches werden als homogen-poröse Körper modelliert. Nicht nur die Porosität, sondern auch die Steifigkeit, Lösbarkeit, Diffusionseigenschaften, der Anteil chemisch instabiler Stoffe und andere Parameter müssen bekannt sein, um die klinisch relevanten Effekte solcher modellierten Implantate oder Deformationen vollständig bestimmen zu können. Entsprechende Stofftransport-Simulationen werden, gemäß dem neuesten Stand der Technik, einschließlich Diplomarbeiten am LSTM-Erlangen, mit Finite-Differenzen-Verfahren auf demselben Rechengitter wie LBM und mit dem LBM-Löser direkt, bei jedem Zeitschritt gekoppelt, mit Verwendung von Simulationsschemata und Randbedingungen zweiter Ordnung (also mit der LBM abgestimmt) durchgeführt.

Für jede dieser Arten von Objekten sollten deshalb schon im voraus Datensätze mit repräsentativen Werten vorhanden sein. Dessen Entwicklung ist kein Bestandteil der vorliegenden Beschreibung. Auch die Auswahl optimaler Modelle für poröse Medien läßt sich nicht eindeutig beschreiben und soll während der Erzeugung der Datensätze bestimmt werden. Sie müssen aber die Struktur der Navier-Stokesschen Gleichungen im wesentlichen behalten, so daß ihre Berechnung mit LBM mittels einer Einführung entsprechend berechneter "Modellkräfte" und "Modellspannungen" möglich bleibt. Die geforderten Angaben von Diffusionseigenschaften sollen eine Berechnung von Fluktuationen in der Temperatur und in der Konzentration von wichtigen Stoffen und Blutbestandteilen erlauben. Entsprechende Diffusions-Konvektions-Gleichungen werden in Zusammenhang mit LBM-Simulationen am besten mittels angepaßter Finite-Differenzen-Verfahren berechnet.

Rechengitter: Auch wenn die ursprüngliche Grauwert-Datei nicht isotrop ist, d.h. wenn der Schritt in der axialen Richtung nicht gleich dem in den Schnittebenen ist, soll die binär segmentierte Geometrie auf einem isotropen Gitter dargestellt werden. Ein Vorteil davon ist, daß numerische Fehler vermieden werden, die aufgrund von Gitteranisotropien bei allen Rechenverfahren auftreten und sich bei geringer Auflösung als wesentlich erweisen. Ein weiterer, eng verbundener Vorteil ist, daß einfachere und genauere Varianten der Lattice-Boltzmann-Methode eingesetzt werden. Die Übertragung auf ein isotropes Gitter kann sehr effizient mit Filtern und Verfeinern verbunden werden. Bei der stark gekrümmten und verzweigten Gefäßgeometrie, die in den meisten Fällen behandelt werden soll, sind diese Vorteile von großer Bedeutung. Wenn das Filtern wesentliche Änderungen der Topologie zur Folge hat, wird dies signalisiert und dem Benutzer (z.B. dem überwachenden Arzt) wird eine einfache Möglichkeit gegeben, die Filtereigenschaften zu modifizieren oder eine der entstandenen Topologien als "richtig" auszuwählen. Das Filtern soll ortsabhängig sein und im Bereich der untersuchten Gefäßdeformation einem (nicht-isotropen, lokal angepaßten) Glättungsoperator entsprechen, üblicherweise mit einer Filterbreite von mindestens 3 und maximal 12 Voxeln des ursprünglichen Gitters, um nicht-physikalische "Wellen" der Gefäßwand im Bereich der Deformation zu vermeiden. Die räumliche Auflösung des entstehenden isotropen Datengitters soll da bei wenigstens 15 bis 20 Voxel entlang der kürzesten Dimension von "normalen" Gefäßen in der Nähe von der untersuchten Gefäßdeformation gewährleisten.

Dieselbe minimale Auflösung gilt auch für kleinere wichtige Gefäße und sogar besonders bei Stenosen, da in dem verengten Lumen höhere Geschwindigkeiten und Gradienten auftreten. Nur bei kleineren Gefäßen mit automatisch bestimmten und typischerweise geringem Ausfluß kann die Auflösung bis auf 6 bis 12 Gitterpunkte (je nach der Variante des verwendeten Lattice-Boltzmann-Schemas) im Durchmesser vermindert werden. Noch kleinere Gefäße (Durchmesser unter 30% von 18 Voxel) werden modelliert und nicht simuliert. Dabei wird eine pulsierende laminare Strömung angenommen und eine dritte Datenbank zum Einsatz gebracht, die eine ausführliche Menge von berechneten kurzstreckigen Einflußströmungen in Gefäßen mit verschiedener Krümmung enthält. Das Geschwindigkeitsprofil am Einströmungsrand kann bei diesen Berechnungen nur mit 2D-Polynomen niedriger Ordnung angegeben werden, da die Aufgabe darin besteht, Geschwindigkeitsdaten mit kleinem Gradienten, die aus Flächen mit Durchmesser unter 5 Voxel stammen, in das modellierte Gefäß "hineinzuextrapolieren".

Eine Gitterzerlegung mit anschließender Eliminierung der Teile, in denen keine durchströmten Voxel vorliegen, und einer Verfeinerung (oder auch Vergröberung, um unnötig hohe Auflösung zu vermeiden) soll nach der Erstellung des isotropen Gitters erfolgen. Jedes Gitterteil soll ein Quader sein, vernetzt von einem eigenen isotropen Gitter, dessen minimale Größe von der technischen Spezifikation des jeweiligen Rechners abhängig ist, aber 12 Voxel in jeder Koordinatenrichtung nicht unterschreiten soll, was mit der Effizienz auf Cache-Architekturen, aber auch mit der Auflösung im Querschnitt der meisten behandelten Gefäßen verbunden ist. Die ursprüngliche Unterteilung ist in etwa gleich großen (plus/minus ein Voxel in jeder Richtung) Blöcken, die sich um ein Voxel (des ursprünglichen isotropen Gitters) überlappen. Nach der Verfeinerung bzw. Vergröberung darf das Verhältnis der Gitterschritte in benachbarten Blöcken nur 1:1 oder 2:1 bzw. 3:1 sein. Dies vereinfacht die Kommunikation und erlaubt eine genaue und effiziente Interpolation zwischen Nachbarteilen. Verfeinerte Blöcke werden automatisch in eine entsprechende Zahl von neuen Blöcken unterteilt, falls die für das vorhandene Rechensystem empfohlene Gittergröße wesentlich überschritten wird. Das so beschriebene Verfeinerungsverfahren wird dann rekursiv wiederholt in Blöcken, wo dies notwendig ist, bis alle Gefäße gut aufgelöst sind. Es werden auch die Beschreibungen aller ursprünglichen Blöcke und aller Zwischenstufen-Blöcke bei der Verfeinerung bzw. Vergröberung gespeichert. Dies ist sowohl bei der Visualisierung und Archivierung der Simulationsergebnisse im Algemeinfall wie auch bei der Simulation selbst, falls LBM-MG verwendet werden, notwendig.

Darstellung der Gefäßwände: Die Identifikation der Lage von Gefäßwänden aus den 3D-Tomographiedaten wurde bereits beschrieben. Diese Lage muß für die DNS in einer in der Handhabung leichten Datenstruktur behalten werden. Die einzige Ausnahme tritt auf, wenn das Bounce-back-Verfahren verwendet wird, als Näherung der Haft-Randbedingungen an Gefäßwänden, die als fest und unbeweglich angenommen werden und nur ungenau bekannt sein dürfen. Dies ist eine schnelle, einfache und robuste Variante des Verfahrens, das auch eine exakte Massenerhaltung garantiert. Es ist aber nur bedingt anwendbar, etwa bei relativ gut aufgelösten Gefäßen, die aber fast unbeweglich sind (wie z.B. im Gehirn) oder bei Gefäßen, dessen Bewegung unwichtig für die Dynamik im Gebiet von Interesse ist. Die Wandinformation wird auf das starre Rechengitter abgebildet und mit diesem gekoppelt. In allen anderen Fällen muß die Wandlage gespeichert und gegebenenfalls auch stets erneut (z.B. bei den starken Pulsationen einer AAA) werden. Wenn die Wanddynamik berechnet werden soll, muß ein mechanisches Model der nichtlinearen Elastizität und (bei großen Abweichungen oder starken Belastungen, wie z.B. in größeren Aneurysmen) der Plastizität angegeben werden. Ein Beispiel solcher Modelle ist in Holzapfel & Ogden (a.a.O., 2003) angegeben.

Es werden zwei Möglichkeiten zur Darstellung der Wandlage berachtet. Bei beiden handelt es sich um eine "genaue" oder "scharfe" Darstellung, was die eindeutige Berechnung der Wandbelastung erlaubt und die effiziente Simulation der Wanddynamik erleichtert, aber auch eine Kommunikation dynamischer Information zwischen dem "Wandgitter" und dem für die Strömungssimulation geeigneten 3D-Rechengitter verursacht. Diese Kommunikation kann mittels des bekannten "immersed-boundary" Verfahrens oder einer seiner neueren Varianten wie des "immersed-interface" Verfahrens (siehe oben) erfolgen.

Die erste Wanddarstellungsmethode ist bei stets kleinen Wandabweichungen einfacher zu behandeln: Die Verbindungen zwischen benachbarten Gitterpunkten im isotropen 3D-Rechengitter werden gekennzeichnet, wenn sie eine Gefäßwand überschneiden, wobei auch die Lage des Schnitts auf dem 3D-Rechengitter gespeichert wird. Fälle, in denen ein Gitterpunkt mit einem Überschneidungspunkt (fast) überlappt, werden durch "generische kleine Deformationen", die lokal angebracht werden können, systematisch vermieden.

Die zweite Darstellungsmethode ist klassisch: Es wird ein unstrukturiertes Gitter (siehe 4) automatisch aufgebaut, z.B. mit dem Marching-Cubes-Verfahren oder anderen Standardverfahren der Bildbearbeitung und Bildgebung. CAD-Darstellungen sind auch möglich, aber in der Regel zu teuer zu erzeugen. Sie können z.B. bei der Generierung von Testgeometrien sinnvoll verwendet werden. Standardformats der größten CAD-Softwarehersteller sowie auch weitere Standardformats wie z.B. STL, werden aus diesem Grund einlesbar und in eine interne Darstellung der Oberflächen-Geometrie automatisch umzuwandeln sein. Das Verhältnis von jedem Punkt dieser Darstellung zum 3D-Rechengitter ist aufgrund der maximal vereinfachten Struktur des isotropen Rechengitters leicht zu ermitteln, auch wenn sich das Wandgitter nach jedem Zeitschritt bewegt hat. Diese Methode ist für größere Wandabweichungen, wie z.B. ein großes Aneurysma oder Teile der Aorta, gut geeignet.

Bei der "virtuellen Operation" von Gefäßdeformationen wird der räumliche Verlauf der Gefäßwände bezüglich des 3D-Rechengitters verändert und gegebenenfalls werden auch Implantate zur Geometriebeschreibung zugefügt, indem ihre "Skelette" (siehe Geometriebestimmung) erst an den Wandverlauf angepaßt und danach aufgebereitet (siehe oben) werden oder indem interaktiv ausgewählte Teilbereiche des durchströmten Gefäßvolumens nicht mehr als Blut, sondern als poröses Medium (s.o.) gekennzeichnet wird. Zur Änderung der Wandgeometrie wird in jedem Fall ein unstrukturiertes 2D-Gitter manipuliert. Gegebenfalls wird das Endergebnis der Geometrieänderung von der Form eines deformierten 2D-Gitters in eine neue Kennzeichnung der wanddurchdringenden Gitterpunktverbindungen auf dem fixierten 3D-Rechengitter umgerechnet. Für die Manipulation des unstrukturierten 2D-Gitters wird ein Teil seiner Punkte interaktiv ausgewählt und danach als elastisch verbunden gekennzeichnet. Dessen Ausgangslage wird als Referenzkonfiguration betrachtet, in der keine unausgeglichenen Spannungen vorliegen. Der Rest der Punkte bleibt fixiert.

Die Möglichkeit, Wandbewegungen nachzubilden und mit pulsierender Strömung nach der Lattice-Boltzmann-Methode (LBM) gekoppelt zu betrachten, stellt einen entscheidenden Vorteil der Erfindung dar. Dabei können die Wandbewegungen synthetisch aus 2D-Darstellungen, insbesondere als Gitter, ebenso aus mehreren Tomographieaufnahmen ermittelt werden. Die modernen Tomographiegeräte lassen eine Aufnahme zeitperiodischer Bewegungen des Herzens oder von größeren Aneurysmen zu. Über die Kopplung mit der Strömung kann die mechanische Belastung der Wand bestimmt werden, ohne daß irgendein invasiver Eingriff vorgenommen werden muß. Diese Vorgehensweise ist genauer und kostengünstiger als direkte Aufnahmemethoden (Ultraschall, Sine-Konstrast-MRI) für die Erfassung der Blutströmung und erlaubt eine zuverlässige Einschätzung der Wandspannungen. Bei der Bearbeitung einer Zeitreihe von Aufnahmen wird ein Zustand, der dem am Ende der Diastole vorgefundenen nahe liegt, als Referenzkonfiguration erzeugt.

Durch grafisch interaktive Bewegung einzelner "elastischer" Punkte in der selbstständigen Wanddarstellung kann die Lage aller anderen Punkte verändert werden. Nach jeder Bewegung werden einander sehr nahe gekommene Punkte in einem Punkt vereinigt und es wird danach die Möglichkeit gegeben, eine automatische Optimierung der Anzahl und Lage der "elastischen" Punkte auf dem 2D-Gitter zu veranlassen. Veränderungen in der Topologie der gesamten Gefäßgeometrie werden anhand der bereits vorhandenen Achsenstruktur nach jeder Bewegung automatisch erkannt und gleich signalisiert.

Verfahren zur "Aufblasung" der Implantate, unter Einbeziehung ihrer elastischen Eigenschaften, bis sie gegen die Gefäßwände fest gedrückt sind, scheinen nicht in geeigneter Form in der Literatur berichtet zu sein. Dessen Prinzip soll die numerische Berechnung der Verdrängung und Deformierung des Implantats aus einem Ausgangszustand sein, die von der Referenzform des Implantats und einer Ausgangslage bestimmt wird, die interaktiv von einer Bedienungsperson, z.B. vom Arzt, vorgegeben wird. Sobald das Implantat nach einer Bewegung in die Richtung einer Gefäßwand sie an einem Punkt erreicht hat, wird es als hyperelastisches bzw. elasto-plastisches Objekt (als 3D-Körper, als Hülle oder als Vernetzung verknüpfter "Drähte") dargestellt. Während der weiteren Anpassung der Lage des Implantats wird die Gefäßwand bei typischen Anwendungen unbeweglich und steif gehalten. Bei der Ausbreitung des "Ballons" wird die benötigte Kraft (gegen den Widerstand des Implantats) ermittelt und es wird eine maximale (verstellbare) Kraft und der minimale Gefäßdurchmesser im Ballonbereich nicht überschritten. Am Ende wird der "virtuelle Ballon" entfernt. Das gesamte Verfahren der Plazierung eines Implantats schließt keine Berücksichtigung des im Gefäß vorhandenen Blutes ein. Die entstandene Geometrie wird weiter genau so behandelt, wie bei bereits tatsächlich im Patienten plazierten Implantaten, deren Lage wie oben beschrieben registriert worden ist.

Strömungssimulation: Die Blutströmung wird als inkompressibel modelliert. Die Bestimmung des Drucks wird dennoch durch zwei Faktoren erschwert: Einerseits erlaubt die Elastizität der Gefäßwände die Ausbreitung von Druckwellen, die in der Standardformulierung inkompressibler Strömungslöser-Algorithmen nicht behandelt werden können. Andererseits stellt jede sehr komplexe Geometrie eine Voraussetzung für Schwierigkeiten bei der Lösung der Poisson-Aufgabe für die augenblickliche Druckverteilung dar. Z.B. ist der Druck auf Gebieten mit "Löchern" nicht eindeutig definiert. Dieses Problem wird mit dem Einsatz von Lattice-Boltzmann-Lösern vermieden, da der Druck dabei mittels numerischer "Druckwellen", also qualitativ gemäß der eigentlichen Physik, und nicht mittels Poisson-Lösern errechnet wird. Vergleiche zwischen LBM und FVM bei der Lösung von laminaren Strömungen in einem "generischen porösen Medium" (siehe oben) haben gezeigt, daß die LBM-Lösung für den Druck wesentlich besser konvergiert. Dabei haben die "inkompressiblen" Varianten der bestens bekannten LBM technische und weitere Vorteile. Das Problem der elastischen Wände wurde bereits diskutiert.

Darüber hinaus kann die nicht-Newtonsche Rheologie des Blutes auch einen maßgeblichen Einfluß haben, besonders in Strömungsgebieten mit sehr niedriger Blutgeschwindigkeit, z.B. innerhalb von großen Aneurysmen oder in Kapillaren. In großen Adern dagegen ist eine Annahme Newtonscher Rheologie nahezu problemlos. Eine Ausnahme sind Fälle, in denen die Wandspannung untersucht werden soll. Ein etabliertes Modell für die Blutrheologie ist bisher nicht bekannt, obwohl mehrere Messungen unternommen worden sind. Besonders bei der Modellierung von teilweise durchströmbaren Teilen des Gefäßlumens als poröse Medien zeichnet sich die Notwendigkeit ab, Tensormodelle für die Viskosität (sowie die Diffusion unterschiedlicher Blutanteile und mitgetragener Stoffe) anzuwenden.

Bisherige LBM-Anwendungen haben stets skalare (auch nicht-Newtonsche, z.B. turbulente) Viskositätsmodelle zum Einsatz gebracht. Dennoch gibt es noch wenig erprobte Ansätze, die den Einsatz von Tensormodellen erlauben, z.B. in dem der Spannungs tensor in die Definition der LBM-Gleichgewichts-Funktion einbezogen wird. Wenn es zur Strömungsnachbildung erforderlich wird, auf komplizierte Spannungsensoren zurückzugreifen, kann es auch notwendig werden, Information aus benachbarten Voxeln zu verwenden und Finite-Differenzen zu berechnen.

Die vollkommene Lokalität bei allen LBM-Algorithmen führt zu einer hervorragenden Skalierung auf Parallelrechnern. Die Anzahl von Operationen und von Freiheitsgraden pro Gitterpunkt ist moderat, so daß auf Prozessoren mit ausreichendem "schnellen Speicher" (Vektor-Register bzw. Cache-Levels) eine hohe Rechengeschwindigkeit erreicht wird. Die Lokalität erlaubt sogar die Berechnung einer guten Näherung des gesamten Spannungstensors nur aus Daten, die auf dem jeweiligen Gitterpunkt lokal vorhanden sind. Besonders in Wandnähe ist dies ein Vorteil vor FVM oder Finite-Differenzen-Methoden (FDM), da die genaen Position und Bewegung der Wand nicht direkt in der Spannungsberechnung benötigt wird. Es liefert auch einen Effizienzvorteil bei der Berechnung von Spannungen. Wenn die FDM-Näherung dennoch berechnet wird, kann sie durch einen Vergleich mit dem Ergebnis der "lokalen" Berechnung der Spannungen einen empfindlichen Test der Genauigkeit in der LBM-Berechnung bei höheren Reynoldszahlen liefern. Der vollexplizite Charakter der LBM bedeutet, daß alle Berechnungen mit dieser Methode einen ausreichend kleinen Zeitschritt erfordern (siehe oben). Andererseits stellt die Notwendigkeit, den Pulsverlauf verläßlich zu diskretisieren, eine ähnliche Anforderung. Der Vergleich beider Einschränkungen ist von dem Erfinder als eine Bestätigung dafür empfunden worden, daß der explizite Charakter der Rechnungen bei Blutgefäßen keine wesentliche Rolle spielt. Die Auflösung, die bei den steilen Zeitgradienten des Druckverlaufs an verschiedenen Orten mit entsprechender Verzögerung erforderlich ist, bedingt die Auswahl des Zeitschritts auf einen konstanten und kleinen Wert. Damit können pulsierende Strömungen, bei denen die Druckschwankungen einen komplexen Zeitverlauf haben, dargestellt werden und insbesondere feine und schnelle Bewegungen im Fluid miterfaßt werden. Implizite Verfahren auf der Basis von LBM (aber meistens mit einer räumlichen Diskretisierung durch FEM oder FVM und nicht, wie bei den Standard-LBM, daß auch in dem erfindungsgemäßen Verfahren vorgesehen ist, durch FDM) sind für statische Probleme entwickelt worden. In dem hier betrachteten Fall haben diese keine Vorteile. Stattdessen wird hier bei stationären Aufgaben auf LBM-MG gesetzt.

Post-Processing: Die wichtigsten Größen bei der Einschätzung der mechanischen Belastung von Gefäßwänden sind die von dem Blutstrom übertragenen Spannungen. Alle Komponenten des Spannungstensors, differenziert in einen viskosen und einen Druckanteil, werden bei DNS mit LBM an jedem Gitterpunkt innerhalb des durchströmten Gebiets zu jedem Zeitschritt als Teil des Simulationsverfahrens und nicht als zusätzliche Rechenlast unabhängig von anderen Punkten und Zeiten berechnet, was die Erstellung von Zeitreihen an voraus von dem Arzt bestimmten Kontrollpunkten erleichtert. Anhand dieser Daten können nach Beendung der Simulation Statistiken wie z.B. ein Oszillationsindex für die Scherspannung (als Vektor tangential zur Wand) sehr schnell errechnet und dargestellt werden. Eine Berechnung von Statistiken während der Simulation verursacht einen wesentlichen Zeitaufwand und soll nur dann eingeschaltet werden, wenn räumlich gemittelte Statistiken, wie z.B. der globale Massenstrom durch einen gewissen Querschnitt oder effektiv die Verweilzeit in einem ausgewählten Volumenteil des durchströmten Gebietes gewünscht werden.

Neben den Spannungen können auch weitere hydrodynamische Felder, wie die 3D-Vektorfelder der Geschwindigkeit oder des Wirbels (berechnet als Anwendung einer Finite-Differenzen-Näherung des Rotationsoperators an das Geschwindigkeitsfeld), wichtige Informationen liefern. Da die Strukturen dieser Felder einem Strömungsmechaniker vertraut und für diesen aussagekräftig sind, einem Arzt aber mit großer Wahrscheinlichkeit von wenig Bedeutung erscheinen werden, müssen ihre Auswirkungen auf die Wandspannungen, Verweilzeit, etc. dem Arzt auf Verlangen mit einfachen und leicht verfügbaren Erklärungen vertraut gemacht werden können. Z.B. sollen starke Wirbel als Gebiete mit niedrigem relativen Druck innerhalb des Lumens und potentiell hohen Wandscherungsraten bezeichnet werden oder sollte die Länge von Stromlinien als Indikator möglicher Verklottung oder Plaquenbildung benutzt werden, usw. Einstellbare Schwellwerte und die automatische numerische Angabe der abschätzbaren Belastungen erleichtern oder gar ermöglichen erst für einen Arzt die Identifizierung und Auswertung der Gefahr individueller Gefäßdeformationen.

Zusätzlich sollen auch die für CFD- und Visualisierungs-Software üblichen Darstellungsmodalitäten verfügbar sein, wie z.B. die mit einem Beispiel vorgestellten 3D-Isoflächen von Geschwindigkeit (oder anderen Vektorfeldern wie dem Druckgradienten und Skalarfeldern wie der Konzentration eines Präparats im Blut), Strömungslinien, Streichlinien und LIC-Verfahren (Linienintegralfaltung). Lokal räumlich und zeitlich (in-phase) gemittelte Daten können leicht auf das ursprüngliche DICOM-Raster projiziert werden. Der pulsierende Anteil der Geschwindigkeit und des Wirbels in ausgewählten Schnitten können ohne großen Aufwand während der Simulation errechnet und danach visualisiert werden.

Darstellungen der Ergebnisse werden immer in einem Standardformat (DICOM) gespeichert, um die automatische Dokumentierung, aber auch die 3D-Visualisierung (durch Rotation, Zoomen, etc.) mit bereits vorhandener und vertrauter Software für die medizinische Bildgebung zu erleichtern. Die Auflösung der Ergebnisdaten muß vor der Simulation gewählt werden und soll nicht gröber als die der des ursprünglichen 3D-Tomographie-Datensatzes sein. Die während der Simulation eigentlich in der Zeit fortbewegten Felder sind viel zu groß und dazu auch nicht direkt in der medizinischen Praxis brauchbar, um ihr Speichern zu rechtfertigen.

Die Erfindung wird im folgenden weiter anhand von Beispielen und der Zeichnung erläutert. In der Zeichnung zeigen:
1 eine Prinzipdarstellung des erfindungsgemäßen Verfahrens,
2 eine schematische Darstellung des ausschließlich auf rechteckigen Gittern und in der Form von Voxelinformation erfolgenden automatisierten 3D-Datenflusses bei dem erfindungsgemäßen Verfahren,
3 ein Beispiel einer binär-segmentierten realen Gehirngefäßgeometrie, bei dem das im Zentrum liegende Aneurysma deutlich erkennbar ist,
4 eine automatische Vernetzung der Oberfläche der in 2 dargestellten Geometrie mit einer unstrukturierten 2D-Datei mittels einen in kommerzieller Software implementierten Marching-Cubes-Algorithmus,
5 eine Momentaufnahme der Größe einer Geschwindigkeitskomponente für ein dargestelltes Gefäß, dargestellt als Isoflächen, die einem positiven ("hellgrau"), einem negativen ("dunkelgrau") und dem Nullwert ("gepunktelt grau") entsprechen.
6 eine Darstellung eines augenblicklichen Strömungszustandes, bei der die getönten Linien tangential entlang dem Geschwindigkeitsvektorfeld verlaufen und die Größe des Geschwindigkeitsvektors lokal durch die Tönung kodiert ist,
7 Staugebiete, Aufwirbelungsgebiete und ein Laminargebiet einer Gefäßstruktur, dargestellt durch ein in kommerzieller Software implementiertes LIC-Verfahren,
8 den Zeitverlauf von Geschwindigkeitskomponenten entlang der drei im Raum und in der Zeit fixierten Koordinatenrichtungen des kartesischen Rechengitters an einem Kontrollpunkt innerhalb eines Aneurysmas und
9 den Zeitverlauf von Scherspannungskomponenten an einem Kontrollpunkt innerhalb eines Aneurysmas wie in 9, wobei die nicht-diagonalen Elemente des hydrodynamischen Spannungstensors dargestellt sind.
1 zeigt eine schematische Darstellung des erfindungsgemäßen Verfahrens. Eine der Hauptvoraussetzungen für die höhere Effizienz des Verfahrens im Vergleich mit bisherigen Ansätzen ist die Beschränkung aller auftretenden 3D-Datenstrukturen auf rechteckige gleichmäßige Gitter, was in 2 schematisch dargestellt ist. Für 1D- (z.B. Raumkurven) und 2D- (z.B. Gefäßwanddarstellungen) Daten gelten keine derartigen Beschränkungen. Dadurch wird die Effizienz kaum beeinträchtigt, da die Daten niedriger Dimension in schnellen Speicherarten auf gegenwärtigen Rechnern gelagert werden können, und deswegen, besonders bei insgesamt größeren Rechengittern, einen unwesentlichen Anteil der Zugriffs- und Rechenzeit in Anspruch nehmen.
Als Anwendungsbeispiel ist ein Gehirnaneurysma veranschaulicht, das eine ausreichend komplexe Geometrie aufweist. Die Ergebnisse der Schritte des erfindungsgemäßem Verfahrens werden hier einheitlich dargestellt.
3 zeigt eine binär-segmentierte Gehinrgefäßgeometrie mit im Zentrum liegenden Aneurysma, das deutlich erkennbar ist. 4 zeigt die Oberfläche der in 3 erscheinenden Gefäßstruktur, in einer Vernetzung mit einer unstrukturierten 2-D Datei.
Bei der in 5 gezeigten Momentaufnahme einer Geschwindigkeitskomponente ist angesichts der Haftbedingungen an den Ge fäßwänden die gesamte Fläche der Wände in gepunkteltem Grau zu sehen. Das Gefäßinnere ist dunkelgrau dargestellt. Unmittelbar vor dem großen Bogen des Hauptgefäßes wird ein Teil des Stroms in das krankhafte Gebilde abgelenkt. Erwartungsgemäß ist dies nur ein sehr kleiner Teil des gesamten Stromvolumens. Wie 6 veranschaulicht, sind innerhalb des Aneurysmas die Geschwindigkeiten sehr niedrig und ein Staugebiet kann in dem gespitzten Teil des Sackes erkannt werden.
7 veranschaulicht eine auf dem LIC-Verfahren beruhende Darstellung auf einer Oberfläche innerhalb der Gefäßstruktur, deren Abstand von der Gefäßwand regelmäßig klein ist. Die Staugebiete (mit erhöhten Druckbelastungen), Aufwirbelungsgebiete (wo eine relativ höhere Scherspannung für besseren Ausbau der Wand sorgt) und das Laminargebiet sind sehr deutlich erkennbar.
In 8 und 9 ist der Zeitverlauf von hydrodynamischen Größen an einem Kontrollpunkt innerhalb des Aneurysmas dargestellt. 8 zeigt Geschwindigkeitskomponenten und 9 Scherspannungskomponenten. Der Pulsverlauf ist klar erkennbar. Wegen der sehr niedrigen Geschwindigkeiten innerhalb des Aneurysmensacks ist eine genaue Periodizität auch nach 10 Pulsperioden noch nicht erreicht. Der Zeichenwechsel einer der Spannungskomponenten, was auch mit einem standardisierten Skalarindex darstellbar ist, weist auf erhöhte mechanische Belastung.
Der größte Aufwand bei der Anwendung der erfindungsgemäßen Verfahrens ist nicht die Segmentierung und (anders als bei manchen der bekannten Ansätze zur numerischen Blutstromsimulation) auch nicht die Vernetzung der Gefäßgeometrie, sondern die rein rechnerische Leistung zur aufgelösten numerischen Simulation. Die entstehenden großen Datenmengen sind nur nach wesentlicher Bearbeitung und Reduktion für einen praktischen Entscheidungsprozeß verwertbar, dann aber läßt die gute Auflösung weitgehende Schlußfolgerungen über den Zusammenhang zwischen den Zustand der Gefäßwände und der hämodynamischen Belastungen zu.
Die Erfindung läßt sich wie folgt zusammenfassen: Dreidimensionale (3-D) Grauwertdaten aus Tomographie-Aufhahmen (MRT oder Röntgen-CT) werden mit Varianten geeigneter bekannter Algorithmen für medizinische und allgemeine Bildgebung binär segmentiert. Entstandene 3-D Gefäßgeometrien werden für eine aufgelöste numerische Simulation ihrer Durchströmung auf kubischen Rechengittern mit Lattice-Boltzmann-Strömungslösern verwendet. Benötigte Randbedingungen werden aus spezifisch entwickelten Datenbanken oder, insofern machbar, aus individuellen nichtinvasiven Messungen extrahiert. Nur wo und wenn die Wanddynamik wichtig ist, wird sie mittels immersed-boundary, immersed-interface oder ähnlichen Verfahren mitsimuliert. Alternativ kann sie aus zeitaufgelösten (4-D) Tomographiedaten mit denselben Segmentierungs-Verfahren erstellt werden, was einen geringeren Rechenaufwand und durch Verzicht auf eine Modellierung der Wandbeschaffenheit auch einen geringeren Rechenfehler verursacht. Im Patienten vorhandene Implantate, Vasoplastik, Plaquen und ähnliches werden identifiziert und in die Simulationsgeometrie eingebracht, wie auch virtuelle (im Patienten noch nicht entstandene) Geometrieänderungen, Implantate und ähnliches mit rechnerisch unaufgelöster Feinstruktur als elastischer, poröser Körper modelliert. Ausschließliche Verwendung finden rechteckige kartesischen Gitter für 3D-Daten und die Lattice-Boltzmann-Methode erfahren sichert Effizienz (wenige Stunden Rechnung, maximal 1 Stunde Arzteinsatz pro Patient), Standardisierung (Gebrauch erfordert keine nicht-medizinische Ausbildung, Eingaben und Ergebnisse in DICOM-Format) und Genauigkeit.
Zwar ist die Erfindung obenstehend im Zusammenhang mit Blutgefäßen beschrieben worden, sie soll jedoch nicht hierauf beschränkt sein, insbesondere auch andere fluidführende Gefäße mitumfassen.

Claims

Verfahren zur Darstellung von dreidimensionalen Strömungsdaten in fluidführenden Gefäßen des menschlichen oder tierischen Körpers unter Verwendung der Lattice-Boltzmann-Methode, bei dem dreidimensionale Gewebestrukturdaten auf regelmäßigen Datengittern (Geometrie-Voxeldaten) des darzustellenden Gefäßbereichs einschließlich dessen Wandung verwendet werden, dieselben Datengitter oder daraus vollautmatisch erzeugte regelmäßige isotrope Rechengitter (Voxel-Gitter) für die Gefäßdarstellung und für die Bestimmung des Strömungsverlaufs ausgewählt werden und der Verlauf der Fluidströmung in zeitlichen Schritten bestimmt wird, dadurch gekennzeichnet, daß ein oder mehrere zu untersuchende Gefäßabschnitte ausgewählt werden, ein Bereich des Gefäßes aus den dreidimensionalen geometrischen Daten ausgewählt wird, der den oder die Gefäßabschnitte enthält und in dem die Strömungszustände dargestellt werden sollen, die dreidimensionalen geometrischen Daten zur Darstellung der Gefäßwände segmentiert werden, die Ein- und Ausströmrandflächen des Gefäßbereiches festgelegt werden, die Gitterrastergröße des Rechengitters abhängig von den Abmessungen des ausgewählten Gefäßbereiches und den Parametern der darin stattfindenden Strömung ausgewählt wird, das Rechengitter zerlegt wird und die Teile des Gitters eliminiert werden, die keine durch- oder angeströmten Voxel enthalten, die zeitliche Bestimmung der Fluidströmung so lange durchgeführt wird, bis die Strömungsdarstellung ausreichend genau ist, die für die Darstellung vorgesehenen Charakteristiken der Strömung und der mechanischen Belastung der Gefäßwände aus der Strömungslösung bestimmt werden und die Ergebnisse in für die Bildgebung oder Archivierung optimierte Datenstrukturen übertragen und gespeichert werden.
Verfahren nach Anspruch 1, dadurch gekennzeichnet, daß die Auflösung des Rechengitters in Gebieten erhöht wird (Verfeinerung), in denen eine größere Zahl von Gitterpunkten benötigt wird und eine erneute Segmentierung durchgeführt wird.
Verfahren nach Anspruch 1, dadurch gekennzeichnet, daß die Auflösung des Rechengitters in Gebieten herabgesetzt wird (Vergröberung), in denen eine geringere Zahl von Gitterpunkten benötigt wird.
Verfahren nach einem der Ansprüche 1 bis 3, dadurch gekennzeichnet, daß ein für dreidimensionale Lattice-Boltzmann-Berechnungen entwickltes Mehrgitter-Verfahren (Multigrid) zur Beschleunigung verwendet wird.
Verfahren nach einem der Ansprüche 1 bis 4, dadurch gekennzeichnet, daß ein blockstrukturiertes Gitter verwendet wird, dessen Blöcke sich überlappen, wobei in benachbarten Blöcken die Gitterschritte im Verhältnis 1:1, 2:1 oder 3:1 sind.
Verfahren nach einem der Ansprüche 2 bis 5, dadurch gekennzeichnet, daß die Verfeinerung rekursiv und bedarfsweise in Blöcken wiederholt durchgeführt wird, bis alle Gefäße gut aufgelöst sind.
Verfahren nach einem der Ansprüche 1 bis 6, dadurch gekennzeichnet, daß ein Gefäßabschnitt ausgewählt wird, dem ein Gefäßstrukturteil zugefügt werden soll.
Verfahren nach Anspruch 7, dadurch gekennzeichnet, daß Angaben zur Geometrie des Gefäßstruk turteils, z.B. aus Tomographie- oder CAD-Daten, als Voxeldaten oder als unabhängige Datenstruktur umgewandelt werden.
Verfahren nach Anspruch 7, dadurch gekennzeichnet, daß die dreidimensionalen geometrischen Daten des Gefäßstrukturteils aus einer gespeicherten Gefäßstrukturtabelle oder -datenbank ausgewählt werden.
Verfahren nach einem der Ansprüche 7 bis 9, dadurch gekennzeichnet, daß die Lage des zu untersuchenden Gefäßstrukturteils bestimmt wird und dessen dreidimensionale geometrische Daten in die Strukturdaten des darzustellenden Gefäßbereichs eingefügt werden.
Verfahren nach einem der Ansprüche 7 bis 10, dadurch gekennzeichnet, daß die dreidimensionalen geometrischen Daten des Gefäßstrukturteils oder eines Teils von diesem als poröse Struktur(en) modelliert worden sind.
Verfahren nach einem der Ansprüche 7 bis 10, dadurch gekennzeichnet, daß der Gefäßstrukturteil ein künstliches, natürliches oder durch Bioengineering erzeugtes Implantat oder eine Vasoplastik ist.
Verfahren nach einem der Ansprüche 7 bis 11, dadurch gekennzeichnet, daß der Gefäßstrukturteil eine Gefäßmalformation (z.B. Erweiterung, Verengung) ist.
Verfahren nach einem der Ansprüche 1 bis 13, dadurch gekennzeichnet, daß für die Einström- und Ausström-Randbedingungen aus einer gespeicherten Tabelle oder Datenbank anhand der gemessenen oder geschätzten lokalen Reynoldszahl und der Gefäßkrümmung ein Geschwindigkeitsprofil bestimmt wird.
Verfahren nach einem der Ansprüche 1 bis 13, dadurch gekennzeichnet, daß für die Einström- und Aus ström-Randbedingungen aus einer gespeicherten Tabelle, insbesondere Datenbank, anhand von anatomischen und physiologischen Angaben, z.B. der Lage des Gefäßes im Körper, Druckverlaufsmessungen, Körperzustandsdaten, ein oder mehrere Pulsmuster bestimmt werden und aus diesen der Pulsverlauf in der Nähe des zu untersuchenden Gefäßstrukturteils bestimmt wird.
Verfahren nach einem der Ansprüche 1 bis 6, dadurch gekennzeichnet, daß für die Darstellung der sich möglicherweise zeitlich ändernden Wandungsgeometrie ein unstrukturiertes zweidimensionales Gitter verwendet wird und vor einer Lageänderung von Punkten des zweidimensionalen Gitters die Lage der Schnittpunkte der wanddurchdringenden Gitterpunktverbindungen im Rechengitter gespeichert wird.
Verfahren nach Anspruch 16, dadurch gekennzeichnet, daß das zweidimensionale Gitter mittels bekannter Algorithmen der Bildgebung und Bildverarbeitung, z.B. des Marching-Cubes-Algorithmus, erstellt wird.
Verfahren nach Anspruch 16, dadurch gekennzeichnet, daß die Wandung des Gefäßbereichs separat von dem Gefäßinneren, insbesondere als zweidimensionales Gitter, insbesondere unstrukturiertes Gitter, dargestellt und mit dem für die Strömungssimulation verwendeten dreidimensionalen Rechengitter durch dynamische Gleichungen gekoppelt wird.
Verfahren nach Anspruch 16 bis 18, dadurch gekennzeichnet, daß das zweidimensionale Gitter und das Rechengitter durch die Immersed-Interface-Methode gekoppelt werden und das zweidimensionale Gitter unabhängig vom Rechengitter bei jedem Zeitschritt geändert wird und das Rechengitter dabei nicht modifiziert wird.
Verfahren nach einem der Ansprüche 16 bis 18, dadurch gekennzeichnet, daß die Wandung des Gefäßbereichs dargestellt wird durch Verbinden der benachbarten Gitterpunkte, die die Wandung durchschneiden (wanddurchdringende Gitterpunktverbindungen), wobei die Lage der Schnittpunkte im Rechengitter gespeichert wird.
Verfahren nach Anspruch 20, dadurch gekenn eichnet, daß, während ein Teil der Punkte des zweidimensionalen Gitters fixiert bleibt, die Verbindung der anderen Punkte als hyperelastisch definiert wird und von diesen Punkten die Lage ausgewählter Punkte verändert wird, womit die Lage der mit diesen Punkten elastisch verbundenen Punkte auch verändert wird.
Verfahren nach einem der Ansprüche 16 bis 19, dadurch gekennzeichnet, daß die Gefäßwandung aus mehreren Tomographieaufnahmen (z.B. aus Röntgen-, MRT-, Doppler-Verfahren) rekonstruiert wird, die jeweils mit dem Rechengitter gekoppelt werden. 23, Verfahren nach Anspruch 1, dadurch gekennzeichnet, daß der Transport von Stoffen und/oder die Diffusion unerschiedlicher Blutanteile durch demselben Rechengitter zugeordnete Skalarwerte in den zeitlichen Schritten der Bestimmung der Fluidströmung gekoppelt mit bestimmt wird.
Verfahren nach der Ansprüche 1 bis 23, dadurch gekennzeichnet, daß die Verteilung der Spannung, der Geschwindigkeit und des Drucks im strömenden Fluid gezeigt wird.
Verfahren nach Anspruch 24, dadurch gekennzeichnet, daß für jeden Punkt des Rechengitters der Spannungstensor, aufgeteilt in einen viskosen und einen Druckanteil, bestimmt wird, ohne numerische Differenzierung zu verwenden.
Verfahren nach Anspruch 24, dadurch gekennzeichnet, daß der Wirbel und die Dissipation im strömenden Fluid gezeigt werden.
Verfahren nach einem der Ansprüche 24 bis 26, dadurch gekennzeichnet, daß der pulsierende Anteil der gezeigten Größen direkt, aber auch durch Berechnung von charakteristischen Zahlen über eine Pulsperiode dargestellt wird.
Verfahren nach den Ansprüchen 24 und 27, dadurch gekennzeichnet, daß die Schubspannungen und der Druck, dessen Mittelwerte und Fluktuationscharakteristiken an den Gefäßwänden gezeigt werden.
Verfahren nach einem der Ansprüche 1 bis 28, dadurch gekennzeichnet, daß vor der Bestimmung der Fluidströmung die Auflösung der Darstellung der Strömung festgelegt wird, die mindestens die Auflösung der dreidimensionalen Strukturdaten ist.
Computerprogrammprodukt zur Umsetzung eines Verfahrens nach einem der Ansprüche 1 bis 29.
Anwendung des Verfahrens nach einem der Ansprüche 1 bis 29 auf Blutgefäße.
Anwendung des Verfahrens nach einem der Ansprüche 1 bis 29, insbesondere 16 bis 23, auf elastische blut- oder fluidführende Gefäße mit bewegten Wandungen.
System zur Darstellung von dreidimensionalen Strömungsdaten in fluidführenden Gefäßen des menschlichen oder tierischen Körpers unter Verwendung der Lattice-Boltzmann-Methode, bei dem dreidimensionale Gewebestrukturdaten auf regelmäßigen Datengittern (Geometrie-Voxeldaten) des darzustellenden Gefäßbereichs einschließlich dessen Wandung verwendet werden und dieselben Datengitter oder daraus vollautmatisch erzeugte regelmäßige isotrope Rechengitter (Voxel-Gitter) für die Gefäßdarstellung und für die Bestimmung des Strömungsverlaufs ausgewählt und der Verlauf der Fluidströmung in zeitlichen Schritten bestimmt wird, insbesondere zur Durchführung des Verfahrens nach einem der Ansprüche 1 bis 29, dadurch gekennzeichnet, daß Mittel zum Auswählen von einem oder mehreren zu untersuchenden Gefäßabschnitten vorgesehen sind, Mittel zum Auswählen eines Bereichs des Gefäßes aus den dreidimensionalen geometrischen Daten vorgesehen sind, der den oder die Gefäßabschnitte enthält und in dem die Strömungszustände dargestellt werden sollen, Mittel zum Segmentieren der dreidimensionalen geometrischen Daten zur Darstellung der Gefäßwände vorgesehen sind, Mittel zum Festlegen der Ein- und Ausströmrandflächen des Gefäßbereiches vorgesehen sind, Mittel zum Auswählen der Gitterrastergröße des Rechengitters abhängig von den Abmessungen des ausgewählten Gefäßbereiches und den Parametern der darin stattfindenden Strömung vorgesehen sind, eine Einrichtung vorgesehen ist, die das Rechengitter zerlegt und die Teile des Gitters eliminiert, die keine durch- oder angeströmten Voxel enthalten, eine Einrichtung vorgesehen ist, die die zeitliche Bestimmung der Fluidströmung so lange durchgeführt, bis die Strömungsdarstellung ausreichend genau ist, eine Einrichtung vorgesehen ist, die die für die Darstellung vorgesehenen Charakteristiken der Strömung und der mechanischen Belastung der Gefäßwände aus der Strömungslösung bestimmt, und eine Einrichtung für die Übertragung und Speicherung der Ergebnisse in für die Bildgebung oder Archivierung optimierte Datenstrukturen vorgesehen ist.