Ein ultraschallbasiertes inverses Messverfahren zur Charakterisierung viskoelastischer Materialparameter von Polymeren

Ein ultraschallbasiertes inverses Messverfahren zur Charakterisierung viskoelastischer Materialparameter von Polymeren Von der Fakultät für Elektrotechnik, Informatik und Mathematik der Universität Paderborn zur Erlangung des akademischen Grades Doktor der Ingenieurwissenschaften (Dr.-Ing.) genehmigte Dissertation von M.Sc. Fabian Bause Erster Gutachter: Prof. Dr.-Ing. Bernd Henning Zweiter Gutachter: Prof. Dr. Andrea Walther Tag der mündlichen Prüfung: 24.03.2016 Paderborn 2016 EIM-E/324 iii Kurzfassung Für die numerische Simulation von Schallausbreitungsphänomenen in, und die Dimensionierung von, akustischen (Mess-)Systemen ist die Verwendung realitätsnaher Materialmodelle und Materialparameter ein zentrales Problem. Akustische Materialparameter lassen sich nicht aus quasistatisch ermittelten Kenngrößen, wie sie in Datenblättern der Hersteller von Halbzeugen angegeben werden, berechnen. Bestenfalls können grobe Abschätzungen getroffen werden, wobei diese gerade bei Polymeren für eine fundierte Betrachtung der Schallausbreitung unzureichend sind. In dieser Arbeit wird ein Messverfahren vorgestellt, welches, für eine gegebene polymere Materialprobe, ein komplex-wertiges und frequenzabhängiges Materialmodell, unter Berücksichtigung von Anisotropie sowie Spannungsrelaxationsund Retardationsprozessen, quantifiziert. Die Materialproben werden als hohlzylindrische Wellenleiter ausgelegt und Ultraschall-Transmissionsmessungen zwischen den parallelen Seiten der Probe durchgeführt. Zur Berücksichtigung der Frequenzabhängigkeit werden fünf verschiedene Schallwandler-Paare mit einer aufsteigenden Mittenfrequenz von 750 kHz bis 2,5 MHz verwendet. Jedes der fünf Messsignale trägt nach Durchlaufen der Materialprobe Informationen über die Materialparameter, die räumliche und spektrale Anregung der Probe sowie die Probengeometrie. In einem inversen Verfahren werden diese Informationen voneinander separiert. Die Lösung des inversen Problems erfolgt deterministisch durch iterativen Vergleich einer Vorwärtssimulation des gesamten Messsystems mit den experimentell bestimmten Messdaten. Bei gegebener Lösung des inversen Problems wird ebenfalls eine Abschätzung der Messunsicherheit eines jeden identifizierten Materialparameters durchgeführt. v Abstract For the numerical simulation of acoustic wave propagation in, and the design of acoustic (measurement) systems, the use of reliable material models and material parameters is a central issue. Acoustic material parameters cannot be evaluated based on quasistatically measured parameters, as are specified in data sheets by the manufacturers. At best, rough estimates can be made, which are insufficient for a thorough consideration of acoustic wave propagation, especially in polymers. In this work, a measurement method is presented which quantifies, for a given polymeric material sample, a complex-valued and frequency-dependent material model, taking anisotropy, stress relaxation and creep retardation processes into account. The material samples are designed as hollow cylindrical waveguides. Ultrasonic transmission measurements are carried out between the parallel faces of the sample. To account for the frequency dependency of the material properties, five different transducer pairs with ascending central frequency of 750 kHz to 2,5 MHz are used. Each of the five received signals contains, after passing through the sample, information on the material parameters, on the spatial and spectral excitation of the sample and on the sample geometry, which are separated from each other in an inverse procedure. The solution of the inverse problem is carried out deterministically by iterative comparison of forward simulations of the entire measurement system with the experimentally determined measurement data. For a given solution of the inverse problem, an estimate of the measurement uncertainty of each identified material parameter is calculated. vii Vorwort Die vorliegende Arbeit ist im Rahmen meiner wissenschaftlichen Tätigkeit im Fachgebiet Elektrische Messtechnik der Universität Paderborn entstanden. Ich möchte mich an dieser Stelle bei allen Personen bedanken, die zum Gelingen dieser Arbeit beigetragen haben. Vorangehend zu nennen ist natürlich Herr Professor Bernd Henning, der sich bereit erklärt hat meine wissenschaftliche Arbeit zu unterstützen. Ich danke Herrn Professor Henning für weitreichende Einblicke in die komplexe Welt des Ultraschalls, ein hervorragendes Arbeitsumfeld sowie für die Erstellung des Erstgutachtens für diese Arbeit. Frau Professorin Andrea Walther möchte ich danken, für ihr offenes Ohr bei zahlreichen Diskussionen rund um die Welt der angewandten Mathematik, zentral natürlich meine vielen Fragen zum Algorithmischen Differenzieren und numerischer Optimierung. Auch der gemeinsame Ausflug in die Intervall-Arithmetik bleibt mir sicherlich in guter Erinnerung. Meinen Dank natürlich auch für die Erstellung des Zweitgutachtens. Danken möchte ich auch allen Mitarbeiterinnen und Mitarbeitern im Fachgebiet Elektrische Messtechnik. Insbesondere zu nennen sind Dr. Jens Rautenberg für seine Unterstützung bei der gemeinsamen Bearbeitung des Forschungsprojekts Bestimmung komplexer akustischer Materialdaten, für dessen Förderung ich der Deutschen Forschungsgemeinschaft danken möchte, für seine offene und ehrliche Art und die Unterstützung meiner Laufbahn seit meinem Bachelor-Studium. Dr. Dietmar Wetzlar für ausgedehnte Diskussionen um das Thema Messunsicherheiten. Ich danke Christian Hoof für die Hilfe der bei der Realisierung der Sende-Elektronik, Sabine Schleghuber für die Hilfe bei der Erstellung von technischen Zeichnungen, Inge Meschede für die ganzen organisatorischen Belange und Friedhelm Rump für seine Unterstützung bei der Realisierung mechanischer Aufbauten. Neben dem Team der Elektrischen Messtechnik danke ich auch allen Studenten, die durch ihre Abschlussarbeiten das Voranschreiten des Forschungsthemas mit begleitet haben. Insbesondere möchte ich hier Nadine Feldmann, Manuel Webersen, Leander Claes, Daniel Weber und Thorsten Meyer meinen Dank aussprechen. Dr. Hauke Gravenkamp danke ich für die gute Zusammenarbeit auf dem Gebiet der Simulation von Wellenausbreitung mit der SBFEM. Dr. Boqiang Huang gilt mein Dank für Einblicke in die äußerst interessante Welt der Wavelets und der empirical mode decomposition. Zuletzt gilt mein aufrichtiger Dank meiner Familie, die mich auf meinem Weg stets unterstützt und mit Rücksicht begleitet hat. Meiner Frau Kerstin danke ich auch für ihr Engagement mir die deutsche Grammatik und Orthographie näher zu bringen. Paderborn, den 01.04.2016 Fabian Bause Inhaltsverzeichnis ix Inhaltsverzeichnis Kurzfassung iii Abstract v Vorwort vii Inhaltsverzeichnis xi Symbolverzeichnis xiii 1 Einleitung 1 2 Grundlagen viskoelastischer Materialien 5 2.1 Grundlegende Merkmale polymerer Werkstoffe ................... 5 2.2 Grundbegriffe der Kontinuumsmechanik . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Grundlagen der linearen Viskoelastizität . . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 Analogiemodelle der Makromechanik . . . . . . . . . . . . . . . . . . . . 11 2.3.2 Zeit-Temperatur-Superpositionsprinzip . . . . . . . . . . . . . . . . . . . 23 3 Stand der Forschung zur Materialparameterbestimmung und akustischer Wellenausbreitung 25 3.1 Materialparameterbestimmung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1.1 Auswahl standardisierter Prüfverfahren . . . . . . . . . . . . . . . . . . . 25 3.1.2 Ultraschallbasierte Prüfmethoden . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Modellierung geführter Wellenausbreitung . . . . . . . . . . . . . . . . . . . . . . 31 3.3 Transiente Finite Elemente Simulation in viskoelastischen Medien . . . . . . . . . 32 3.3.1 Grundgleichungen und das Rayleigh-Dämpfungsmodell in der Finite ElementeMethode................................ 33 3.3.2 Zeit-kausale viskoelastische Dämpfung in der Finite Elemente Methode . . 35 4 Wellenausbreitung in viskoelastischen Medien 37 4.1 Implikationen viskoelastischer Medien auf die Wellengleichungen . . . . . . . . . 37 4.2 Viskoelastische Materialmodelle für die Simulation akustischer Wellen in Polymeren 39 4.2.1 Viskoelastische Re-Interpretation der Rayleigh-Dämpfung . . . . . . . . . 39 4.2.2 Diskussion des fraktionalen Zener-Modells . . . . . . . . . . . . . . . . . 43 4.2.3 Generalisierung auf dreidimensionale Betrachtungen und Diskussion der Anisotropie.................................. 45 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter . . . . . . . . . . . 53 4.3.1 Halb-analytischer Ansatz in der Scaled Boundary FEM ........... 54 4.3.2 Methode der Modalen Expansion . . . . . . . . . . . . . . . . . . . . . . 58 xvi Symbolverzeichnis Formelzeichen Einheit Bedeutung EV/m Vektor der elektrischen Feldstärke E0,E1,E2Koeffizientenmatrizen der SBFEM fVN/m3Volumenkraftvektor FN Kraft fHz Frequenz fsHz Abtastfrequenz fcHz Mittenfrequenz G(t)Pa Relaxationsmodul G+(iω)Pa Komplexes Modul G0Pa Speichermodul G0 FZ,∞Pa Grenzwert des Speichermoduls für ω→ ∞ G00 Pa Verlustmodul G00 FZ,∞Pa Grenzwert des Verlustmoduls für ω→ ∞ h(t)Impulsantwort eines Systems H(t)Sprungfunktion (Heavisidefunktion) HWL(iω) Frequenzgang eines Wellenleiters gegeben räumliche Anregung und Empfangscharakteristik H(z)Übertragungsfunktion in der z-Transformierten hV/m=N/C Matrix piezoelektrischer Moduln IEinheitsmatrix IA Elektrische Stromstärke J(t)1/Pa Kriechnachgiebigkeit J+(iω)1/Pa Komplexe Kriechnachgiebigkeit J01/Pa Speichernachgiebigkeit J00 1/Pa Verlustnachgiebigkeit kPMechanischer Kopplungsfaktor (planar) kTMechanischer Kopplungsfaktor (thickness) kλ1/K Thermischer Ausdehnungskoeffizient ktrans Kopplungsfaktor eines Übertragers krad/m Wellenzahl KPa Kompressionsmodul KuSteifigkeitsmatrix (für die FEM) lm Länge eines Körpers LH Elektrische Induktivität mkg Masse mnn-tes gewöhnliches Moment MuMassenmatrix (für die FEM) M0Massenmatrix (für die SBFEM) MN m (Dreh-)Momentenmatrix MS(iω)N/(V m2) Sendefrequenzgang eines Schallwandlers ME(iω)V m2/N Empfangsfrequenzgang eines Schallwandlers NbFormfunktion in der FEM Ni(ζ)Formfunktion in der SBFEM NRMSD % Normalized Root-Mean-Square Deviation p(r,t)Pa Schallwechseldruck Symbolverzeichnis xvii Formelzeichen Einheit Bedeutung pVektor der Ursachen/Optimierungsvariablen ˆqn Amplituden der Knotenpunktkräfte in der SBFEM QmMechanische Güte rm Ortsvektor rim Innenradius RΩOhm’scher Widerstand S1/Pa Matrix der elastischen Nachgiebigkeit (Voigt’sche Notation) SNR dB Signal-zu-Störabstand sσStichprobenstreuung Sxy(iω)Kreuzleistungsdichtespektrum Sxx(iω)Autoleistungsdichtespektrum ts Zeit t0s Signallaufzeit zum L-Puls ∆ts Signallaufzeit zwischen L-Puls und T-Puls tpc m Dicke der Piezoelektrischen Scheibe tem Dicke der Elektrodenschicht auf der Piezokeramik tjs Jitter TK Temperatur TgK Glasübergangstemperatur Tps Pulsdauer Ttot s Totzeit Ts(·)Synchrosqueezed Wavelet Transform UEJ Innere Energie UV Elektrische Spannung ucKombinierte Standardunsicherheit uTypA Typ A Standardunsicherheit uTypB Tpy B Standardunsicherheit UX% Erweitere Unsicherheit bezüglich X% - Vertrauensintervall Um Knotenpunkt-Verschiebungsvektor Vm/s (Schall-)Schnelle Wy(a,b)Wavelet-Transformierte des Signals y Wyy(t,f)Wigner-Ville-Verteilung yVektor der Beobachtungen yp(t)Signalprototyp Y(f,t)Kurzzeitfouriertransfomierte des Signals y(t) ZmRayl m2Akustische Strahlungsimpedanz Zak Rayl Akustische Impedanz ZWΩElektrischer Wellenwiderstand Zin ΩElektrische Eingangsimpedanz xviii Symbolverzeichnis Griechische Buchstaben Formelzeichen Einheit Bedeutung αKs Rayleigh-Dämpfungskoeffizient (∼Steifigkeit) αM1/s Rayleigh-Dämpfungskoeffizient (∼Masse) α(ω)Np/m Schallabsorption (längenbezogen) ˜α(ω)Np/s Schallabsorption (zeitbezogen) α0,L Np/m Longitudinalwellenabsorption bei f0 α0,T Np/m Transversalwellenabsorption bei f0 αth W/(m2K) Wärmeübergangskoeffizient β(ω)1/m Phase der Wellenausbreitung βSV m s/A Kehrwert der Permittivität Γ(·)Gammafunktion γ1/m Wellenausbreitungskonstante δ(t)Delta-Funktion δii Kronecker delta tan δVerlustfaktor εS 33 A s/(V m) Geklemmte Permittivität in 33-Richtung εT 33 A s/(V m) Freie Permittivität in 33-Richtung εMechanischer Verzerrungsvektor [ε]Causchy’scher Verzerrungstensor Kostenfunktion/Zielfunktion ζiLehr’sches Dämpfungsmaß ζLokale Knotenpunktkoordinate ηPa s Dynamische Viskosität θrad Azimutale Koordinate κNp/λAbsorptionskoeffizient pro Wellenlänge λrRegularisierungsparameter λim2/N Eigenwerte der Nachgiebigkeitsmatrix ˜ λ Eigenwert der charakteristischen Funktion der SBFEM µPa Schubmodul für isotropes Material µLPa Schubmodul in Vorzugsrichtung µTPa Schubmodul in transversaler Ebene νFraktionale Ableitungsordnung νG Fraktionale Ableitungsordnung bezogen auf Verluste der Scherung νK Fraktionale Ableitungsordnung bezogen auf Verluste der Dilatation νiso Querdehnzahl /Poisson-Zahl für isotropes Material νLQuerdehnzahl /Poisson-Zahl in Vorzugsrichtung νT Querdehnzahl /Poisson-Zahl in transversaler Ebene ξm Verschiebungsvektor ξnm Knotenpunktverschiebungsvektor (SBFEM) ˆ ξnm Amplitudenvektor der Knotenpunktverschiebungen (SBFEM) Symbolverzeichnis xix Formelzeichen Einheit Bedeutung Φx,y Kreuzkorrelationsfunktion $Gewichtungsvektor ρkg/m3Dichte σN/m2Mechanischer Spannungsvektor σN/m2Mechanischer Eigen-Spannungsvektor [σ]N/m2Causchy’scher Spannungstensor σ(j) zz (r,ω)N/m2Radialer Normal-Spannungsverlauf der Mode j σ(j) rz (r,ω)N/m2Radialer Tangential-Spannungsverlauf der Mode j σtrans Streufaktor eines Übertragers ΣKovarianzmatrix τεs Retardationszeitkonstante τσs Relaxationszeitkonstante τgd s Gruppenverzögerungszeit τs Verschiebezeit τs(a,b)sCandidate group delay ϕm2Skalares Potential Ψm2Vektorpotential Ψ(·)Wavelet ωrad/sKreisfrequenz ωs(a,b)Hz Candidate instantaneous frequency ωe◦Eigen-Winkel 1 Einleitung 1 1 Einleitung Die Entwicklung neuer Werkstoffe wird zunehmend geprägt durch Materialien auf Polymer-Basis, deren Energie-, Ressourcenund Kosteneffizienz stetig zunimmt. Treibende Faktoren dieser rasanten Entwicklung polymerer Werkstoffe sind unter anderem ökologische Zielsetzungen hinsichtlich der Emissionsziele. Beispielsweise ist in den Bereichen Mobilität und Energieerzeugung der Einsatz von Leichtbaustrukturen ein aktuelles und dringliches Forschungsund Entwicklungsthema. Bereits Ende der 1980er Jahre übertraf das (geometrische) Volumen der weltweiten Kunststoffproduktion das Produktionsniveau von Rohstahl [GS11]. Des Weiteren ist die Preisentwicklung verschiedener Metalle ein essenzieller wirtschaftlicher Faktor bei der Entwicklung einer Vielzahl von Produkten. So haben sich die Rohstoffpreise von Kupfer und Zink zwischen 1999 und 2006 mehr als verdreifacht, der Aluminiumpreis stieg in dieser Zeitspanne etwa um den Faktor 1,6 und Eisenerz (Feinerz) hat seinen Preis deutlich mehr als verdoppelt [BGR06]. Hingegen ist die Preisentwicklung polymerer Werkstoffe stabiler, wenn auch eine starke Korrelation zum Rohölpreis vorhanden ist. In Abb. 1.1 ist die Preisentwicklung ausgewählter Thermoplaste und des Rohölpreises über die letzten 10 Jahre exemplarisch dargestellt. Abbildung 1.1: Preisentwicklung exemplarischer Polymere in den letzten 10 Jahren (Quelle: www.kiweb.de). Der zunehmende Einsatz polymerer Werkstoffe führt im Bereich der ultraschallbasierten Messund Prüftechnik sowohl zu neuen Möglichkeiten als auch zu neuen Herausforderungen. So bieten 2 1 Einleitung Kunststoffe erheblich mehr Freiheitsgrade bei der synergetischen Gestaltung konstruktiver sowie akustisch funktionaler Elemente [Rau12]. Nachteilig sind ihre deutlich komplexeren akustischen Eigenschaften. Die frequenzund zumeist richtungsabhängige(n) Schallgeschwindigkeit(en), Absorption 1 und Schallkennimpedanz variieren, anders als bei Metallen, deutlich mit der Temperatur, dem Alter oder dem Wassergehalt. Zudem sind die Absorptionseigenschaften vieler Kunststoffe deutlich komplexer und hinsichtlich der Schallausbreitung erheblich relevanter als bei den meisten Metallen. Ebenfalls spielt der Herstellungsprozess selbst, d.h. das Fertigungsverfahren, die Prozessparameter sowie die Umgebungseinflüsse, eine zentrale Rolle bei der Ausprägung der Materialeigenschaften von Polymeren. Die Bandbreite dieser Schwankungen gilt es im Vorfeld einer jeden Produktentwicklung zu erfassen und bspw. bei der Dimensionierung eines ultraschallbasierten Sensors oder der Beurteilung eines polymeren Prüfkörpers zu berücksichtigen [RBH13]. Einordnung der Arbeit Die quantitative Beschreibung physikalischer Materialeigenschaften erfolgt über Materialmodelle. Im Gegensatz zu Stoffgesetzen ist die Nutzung von Materialmodellen nicht eindeutig, so können für ein Material verschiedene Beschreibungsmodelle herangezogen werden, wobei jeweils deren Gültigkeitsgrenzen zu diskutieren sind. Einem Materialmodell sind zur quantitativen Abbildung von physikalischen Einflussgrößen (Modell-Eingangsgrößen) auf dessen Reaktionsgrößen (ModellAusgangsgrößen) Materialparameter zugeordnet. Für die realitätsnahe numerische Simulation von Schallausbreitungsphänomenen in Messsystemen, d.h. sowohl im akustischen Schallwandler als auch im angrenzenden soliden oder fluiden Messmedium, ist die Verwendung realitätsnaher Materialmodelle und Materialparameter ein zentrales Problem. Ausgangspunkt der vorliegenden Arbeit bildet ein vorangegangenes Forschungsprojekt der Fachgruppe Elektrische Messtechnik der Universität Paderborn zur Charakterisierung stark absorbierender Kunststoffe, welches in der Dissertationsschrift von Dr. Jens Rautenberg mündete [Rau12]. Hier wird ein wellenleiterbasiertes inverses Messverfahren beschrieben, welches mittels einer einzelnen Transmissionsmessung durch eine hohlzylindrische Probe ein reellwertiges transversal isotropes Materialmodell identifiziert. Für die Beschreibung der Schallabsorptionsvorgänge wird das Rayleigh-Dämpfungsmodell genutzt, welches bei transienten Simulationen zwangsläufig nicht-kausales Verhalten zeigt, jedoch noch weit in Wissenschaft und Industrie verbreitet ist. Zudem ist das Modell nur in einem limitierten Frequenzbereich gültig. Die an individuellen Probekörpern identifizierten Materialparameter können dann für die numerische Simulation mit der Finite Element 1 Die Begriffe Dämpfung und Absorption werden im Folgenden unterschieden. Im Kontext der Wellenausbreitung ist der Begriffder Dämpfung die Summe aus Absorption, d.h. Umwandlung in Wärme, und Verlusten durch geometrische Aufweitung der Schallkeule. Im Kontext der Materialmodelle wird zur Modellierung von Verlusten im Sinne von Absorption gebräuchlich ebenfalls von Dämpfung gesprochen. Es wird daher im Kontext der Materialmodelle der Begriffder Dämpfung weiterhin genutzt. 3 Methode (FEM) genutzt werden [Rau12; BGR15b]. Aufgrund der regen Forschungsaktivitäten auf dem Gebiet der Implementierung zeit-kausaler Dämpfungsmodelle für die Elastodynamik innerhalb der FEM, vgl. Kapitel 3, ist es folglich nur konsequent, messtechnische Methoden zur Identifikation der dort genutzten Materialmodelle zu entwickeln, um das Potenzial der neuen Dämpfungsmodelle in Zukunft konstruktiv nutzen zu können. Akustische Materialkenngrößen lassen sich nicht aus quasistatisch ermittelten Materialkenngrößen berechnen. Bestenfalls können hier grobe Abschätzungen getroffen werden. Bezüglich der Symmetrieeigenschaften der Materialmodelle können viele Metalle mit ausreichender Genauigkeit durch ein isotropes Materialmodell (absolute Symmetrie des Materialtensors) beschrieben werden. Die meist geringe frequenzabhängige Absorption kann oft, zumindest bei Betrachtung eines eingeschränkten Frequenzbereichs, mit dem Rayleigh-Dämpfungsmodell hinreichend abgebildet werden, obwohl auch hier bei transienten Simulationen Kausalitätsverletzungen auftreten. Die Materialsymmetrie wird maßgeblich vom Produktionsprozess, im Wesentlichen Extrudieren oder Spritzgießen [CLY92; War75], und ggf. durch die Zugabe von Zuschlagsstoffen wie Kurzglasfasern oder Glaskugeln beeinflusst. Zudem verhalten sich Kunststoffe im Allgemeinen viskoelastisch, d.h. es treten neben der reinen Elastizität der Körper auch Relaxationsmechanismen auf. Ein viskoelastisches Material hat entsprechend eine Art Erinnerungsvermögen an vergangene Belastungen. Viskoelastisches Materialverhalten kann z.B. durch frequenzabhängige komplexe Moduln abgebildet werden. Aus den Realteilen lassen sich die Phasengeschwindigkeiten einer akustischen Welle ableiten, wohingegen die Imaginärteile die Absorption repräsentieren. Der frequenzabhängige Zusammenhang zwischen Realund Imaginärteil der Moduln muss dabei aus Gründen der Kausalität den Kramers-Kronig-Beziehungen [WMM05] genügen. Applikationsbeispiele, in denen die Auslegung des Sensorsystems in hohem Maße von den Materialeigenschaften der eingesetzten Kunststoffe abhängig sind, reichen von der ultraschallbasierten Durchflussmessung, bei der tendenziell Metalldurch Kunststoffrohre substituiert werden, bis hin zur zerstörungsfreien Werkstoffprüfung und Komponentenüberwachung polymerer bzw. hybrider Strukturen, bei denen Kenntnisse über das akustische Verhalten der Polymere als Funktion der Temperatur, des Wassergehalts und der Alterung wichtig sind, um Schadensbilder von natürlichen Veränderungen des Prüfobjekts zu unterscheiden. Weitere Einsatzmöglichkeiten für ein solches Messsystem können sein: •Die Charakterisierung des akustischen Alterungsverhaltens von Kunststoffen (für die Auslegung langlebiger Produkte und Messzellen sowie für die Komponentenüberwachung) • Die gezielte akustische Charakterisierung von Polymerblends (bspw. für die Suche nach geeigneten Kunststoffen zur Impedanzanpassung oder Schalldämpfung) im Kontext akustischer Messund Prüftechnik 4 1 Einleitung • Die Beobachtung der Stabilität der Kunststoffeigenschaften zur Qualitätskontrolle (bspw. Schwankungen der Eigenschaften bezüglich verschiedener Chargen) Wissenschaftliche Fragestellungen In dieser Arbeit wird die Realisierung eines ultraschallbasierten Messsystems zur Charakterisierung viskoelastischer polymerer Werkstoffe beschrieben. Der Hauptanspruch liegt dabei in der Berücksichtigung zeit-kausaler Modellierung transienter Schallausbreitung auf Ebene der Materialgleichung im Frequenzbereich. So ein Materialansatz erlaubt die konsistente Beschreibung von viskoelastischen Polymeren über einen großen Frequenzbereich. Aufbauend auf die am Fachgebiet durchgeführten vorangegangenen Arbeiten (vgl. [Rau12]), soll auch hier das Materialmodell durch die Ausnutzung der Überlagerung von multimodaler geometrischerund Materialdispersion im Wellenleiter in einem inversen Ansatz identifiziert werden. Aus dem Vorhaben lassen sich folgende zentrale wissenschaftliche Fragestellungen formulieren: • Wie lässt sich die Ausbreitung akustischer Wellen in viskoelastischen Polymeren unter der Berücksichtigung der Kausalität beschreiben? Ein Schwerpunkt soll hier auf die konsistente Formulierung eines viskoelastischen Materialmodells unter Annahme einer hexagonalen Materialsymmetrie gelegt werden. • Wie lässt sich so ein Materialmodell in ein effizientes Vorwärtsmodell zur transienten Simulation geführter Wellenausbreitung einbetten? • Wie sind die Schallwandler des Messaufbaus zu entwerfen und zu identifizieren, um im Vorwärtsmodell möglichst realitätsnah abgebildet zu werden? • Wie ist die Sensitivität des Messverfahrens auf die zu bestimmenden Materialparameter zu beurteilen? • Wie ist das nichtlineare inverse Problem zu lösen? Wie können geeignete Startwerte aus den Messsignalen gewonnen werden? • Wie ist die Messunsicherheit der bestimmten Materialparameter schlussendlich zu bewerten? 2 Grundlagen viskoelastischer Materialien 5 2 Grundlagen viskoelastischer Materialien 2.1 Grundlegende Merkmale polymerer Werkstoffe Polymere bestehen aus Makromolekülen, welche wiederum durch Moleküle mit einer großen Anzahl von chemisch gleichen Grundbausteinen charakterisiert sind. Polymere lassen sich grundlegend in Thermoplaste, Elastomere und Duroplaste einteilen [DEE08]. Hauptunterscheidungsmerkmal ist die Art der Vernetzung der Makromoleküle. Während bei Elastomeren und Duroplasten chemische Bindungen zwischen den Makromolekülen entstehen (Vernetzung), basiert der Zusammenhalt bei Thermoplasten auf Verschlaufungen (statistisches Knäuel) und sekundären Bindungen zwischen den unvernetzten Makromolekülen [DEE08]. Daraus ergibt sich aus anwendungsorientierter Sicht ein wichtiger Vorteil von Thermoplasten. Sie sind beliebig oft erweichbzw. schmelzbar, während Duroplaste und Elastomere nach einmaliger Vernetzung nicht mehr löslich sind. Elastomere und Duroplaste unterscheiden sich dadurch, dass Elastomere nur ein schwaches Netz, Duroplaste hingegen ein sehr engmaschiges Netz ausbilden. Der Vernetzungsgrad der Makromoleküle bestimmt letztendlich die temperaturund zeitabhängigen Medieneigenschaften, wie sie im Folgenden betrachtet werden. Je höher der Vernetzungsgrad, um so weniger ausgeprägt ist das Retardationsvermögen des Materials unter Last und um so größer ist die thermische und mechanische Belastbarkeit [DEE08]. Schwerpunkt dieser Arbeit bildet die Untersuchung thermoplastischer Materialen. Thermoplaste können in amorpher (bspw. PMMA, PUR, PS) 2 und in teilkristalliner (bspw. PE, PP, PA, PTFE, PEEK) Form existieren. Bei amorphen Strukturen sind die Makromoleküle statistisch so angeordnet, dass sie einen entropisch günstigen Zustand einnehmen. Teilkristallinität entsteht, wenn sich lineare Makromoleküle regelmäßig anordnen, was bei chemisch gleichmäßigem Aufbau der Ketten vorkommen kann [DEE08]. Die Unterscheidung dieser beiden Gruppen ist gerade hinsichtlich ihres Temperaturverhaltens (und damit ihres Frequenzverhaltens, siehe Kapitel 2.3.2) von großer Bedeutung. Unterhalb der Glasübergangstemperatur Tg eines Polymers befindet sich der Werkstoffim Glaszustand, welcher durch energieelastisches Verhalten charakterisiert ist. Die Glasübergangstemperatur definiert den Übergangsbereich hin zum gummi-elastischen Bereich, welcher 2 Es werden im Folgenden zur besseren Lesbarkeit direkt die Kurzzeichen der Polymere verwendet, vgl. Tabelle der verwendeten Abkürzungen. 12 2 Grundlagen viskoelastischer Materialien basierend auf Federn und Dämpfern charakterisiert werden können [Fer80; Mai10]. Die Struktur der Analogiemodelle definiert das Kriechund Relaxationsvermögen des betrachteten Materials. Aufgrund der Vielzahl verschiedener Modelle, welche zu bestimmten Problemen und Materialen passend erscheinen, sei an dieser Stelle für einen Überblick auf umfangreiche Literatur verwiesen [Mai10; SW00; HN11; CH03; NH13]. Im Folgenden wird nur auf einige exemplarische Modelle eingegangen, welche für die vorliegende Arbeit von Bedeutung sind. ηE (d) η η (a) (b) (c) E E Abbildung 2.4: Elementare Analogiemodelle der Makromechnik: (a) Hooke-Modell, (b) Newton-Modell, (c) Maxwell-Modell, (d) Kelvin-Voigt-Modell. Für einführende Überlegungen ist es sinnvoll, zunächst das Hooke-Modell (Abb. 2.4 (a)) sowie das Newton-Model (Abb. 2.4 (b)) zu betrachten. Die konstitutive Gleichung des Hooke-Modells ergibt sich zu σ(t) = Eε(t);E∈R+. (2.22) Damit ergeben sich die Materialfunktionen durch Betrachtung eines Belastungssprungs in Spannung oder Verzerrung, vgl. Gl. (2.13), zu J(t) = E−1H(t) und G(t) = EH(t) . Wie zu erwarten, werden durch dieses Modell weder Relaxationsnoch Kriechprozesse berücksichtigt. Diese Materialbeschreibung ist demnach nur für ideal elastische Festkörper anwendbar. Über das Boltzmann’sche Superpositionsprinzip, vgl. Gl. (2.15), gelangt man wieder zum allgemeinen Materialverhalten nach dem Hooke’schen Gesetz, in diesem Spezialfall ist dies identisch mit der Materialgleichung des Hooke-Modells. Betrachtet man hingegen die Differentialgleichung (DGL) des Newton-Modells σ(t) = η˙ε(t);η∈R+, (2.23) ergeben sich die Materialfunktionen zu J(t) = tη−1(2.24) G(t) = ηδ(t). 2.3 Grundlagen der linearen Viskoelastizität 13 Das Newton-Modell beschreibt demnach ein ideal viskoses Fluid (Newton’sches Fluid) mit linearem Kriechverhalten und einer unmittelbaren Spannungsrelaxation. Viskoelastische Materialien vereinigen die elastischen und viskosen Eigenschaften. In welcher Form die grundlegenden Federund Dämpfer-Elemente kombiniert werden, hängt dabei von den zu modellierenden Effekten ab. Beispielsweise werden beim Maxwell-Modell (Abb. 2.4 (c)) ausschließlich (exponentielle) Spannungsrelaxationsprozesse aber keine Kriecheffekte modelliert [Mai10]. Das Kelvin-Voigt-Modell hingegen (Abb. 2.4 (d)) bildet nur (exponentielle) Kriecheffekte ab, jedoch keine Spannungsrelaxation. Durch Kombination der genannten Modelle ist es möglich, Materialgleichungen zu erzeugen, welche gleichermaßen Retardation und Relaxation berücksichtigen. (a) (b) E1 E2 η η2 E η1 Abbildung 2.5: Erweiterte Analogiemodelle der Makromechnik: (a) Anti-Zener-Modell, (b) Zener-Modell. Im Folgenden werden zum einen das Anti-Zener-Modell (Abb. 2.5 (a)) und zum anderen das ZenerModell (Abb. 2.5 (b)) betrachtet. Die DGL des Anti-Zener-Modells, mit den Größen η1∈R+ , η2∈R+und E∈R+, lässt sich bestimmen zu σ(t) + η1 E˙σ(t) = (η1+η2)˙ε(t) + η1η2 E¨ε(t). (2.25) Gebräuchlicher ist die Darstellung in der Form σ(t) + τσ˙σ(t) = ˜η(˙ε(t) + τε¨ε(t))(2.26) mit den Substitutionen ˜η=η1+η2;τε=η1η2 E(η1+η2);τσ=η1 E. (2.27) 14 2 Grundlagen viskoelastischer Materialien Damit ergeben sich die Materialfunktionen zu [Mai10] J(t) = J0t+J11−e−t τε(2.28) G(t) = G0δ(t) + G1e−t τσ, mit J0=˜η−1;J1=τσ−τε ˜η(2.29) G0=τε˜η τσ ;G1=˜η(τσ−τε) τ2 σ . An dieser Stelle zeigt sich, dass die Größe τσ die (Spannungs-)Relaxationskonstante und τε die Retardationskonstante beschreibt. Aus physikalischer Perspektive muss der Relaxationsmodul G(t) eine monoton fallende und die Kriechnachgiebigkeit J(t) eine monoton steigende Funktion darstellen. Damit muss für dieses Modell τσ> τε gelten. Abb. 2.6 zeigt schematisch das mit dem Anti-Zener-Modell abbildbare Kriechund Relaxationsverhalten. Wichtige Eigenschaften des Modells sind das unendliche Kriechverhalten mit bleibender Deformation sowie die vollständige reversible Spannungsrelaxation bei instantanem impulsartigen Spannungsanstieg. Die Eigenschaften eignen sich entsprechend eher für die Beschreibung eines viskoelastischen Fluids als für einen viskoelastischen Festkörper. Das Anti-Zener-Modell ist in der Rheologie auch als Jeffreys fluid bekannt [RRP13]. Hingegen ist das Zener-Modell auch als standard linear solid bekannt. Seine Materialdifferentialgleichung mit den Größen η∈R+ , E1∈R+ und E2∈R+ , ergibt sich zu [Zen48; Mai10] σ(t) + η E2 ˙σ(t) = E1ε(t) + E1+E2 E2 η˙ε(t). (2.30) Durch Substitution von E=E1;τε=E1+E2 E1E2η;τσ=η E2, (2.31) ergibt sich die gebräuchlichere Form der DGL σ(t) + τσ˙σ(t) = E(ε(t) + τε˙ε(t)). (2.32) Die Materialfunktionen des Zener-Modells lassen sich dann bestimmen zu [Mai10] J(t) = τσ Eτε +τε−τσ Eτε1−e−t τε(2.33) G(t) = E+Eτε−τσ τσ e−t τσ. 2.3 Grundlagen der linearen Viskoelastizität 15 Wie schon bei der Betrachtung des Anti-Zener-Modells beschreiben τσ die (Spannungs-) Relaxationskonstante und τε die Retardationskonstante mit dem Unterschied, dass mit identischer Begründung wie zuvor nun τσ< τε gelten muss. Die Materialfunktionen des Zener-Modells weisen eine instantane Nachgiebigkeit sowie eine instantane endliche Elastizität auf. Des Weiteren erreicht der Kriechprozess des Materials einen statischen Endwert. Das Materialverhalten ist schematisch in Abb. 2.7 dargestellt. Thermoelastische Interpretation und Generalisierung des Zener-Modells Zur physikalischen Interpretation des Zener-Modells wird die thermoelastische Kopplung betrachtet [Zen48] ε(σ,∆T) = ∂ε ∂σ!T=konst. σ+ ∂ε ∂T!σ=konst. ∆T(2.34) =E−1σ+kλ∆T, welche die Verzerrung eines Körpers neben den anliegenden mechanischen Spannungen auch in Relation zur Temperaturänderung ∆T des Körper setzt, wobei kλ∈R den thermischen Ausdehnungskoeffizienten beschreibt. In dieser Beziehung wird E als isothermales Elastizitätsmodul bezeichnet. Die Änderung der Temperatur sei nun verursacht durch die Verformung des Körpers sowie der Änderung der inneren Energie UE. ∆T(ε,UE) = ∂T ∂ε !UE=konst. ε+ ∂T ∂UE!ε=konst. UE(2.35) Durch die Verformung des Körpers treten zwei Mechanismen auf, welche die Temperaturänderung bestimmen. Durch zeitliche Ableitung von Gl. (2.35) lassen sich diese als Dehnung unter adiabatischen Verhältnissen und Wärmeleitung identifizieren. d dt∆T= ∂T ∂ε !UE=konst. ˙ε+ ∂T ∂UE!ε=konst. ˙ UE(2.36) = d dt∆T!diffus. + d dt∆T!adiab. =−γd˙ε−τ−1∆T Dabei sei τ=τσ in Anlehnung an die Konvention, dass τσ die Relaxation unter konstanter Verzerrung ε beschreibt. Werden die beiden Gleichungen derart kombiniert, dass ∆T entfällt, ergibt sich σ+τσ˙σ=E(ε+τσ(1+kλγd)˙ε), (2.37) 16 2 Grundlagen viskoelastischer Materialien Beliebige Zeiteinheit Spannung bzw. Verzerrung 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 (a) Kriechen (blau) durch sprunghafte mechanische Spannung (rot). Beliebige Zeiteinheit Spannung bzw. Verzerrung 0 50 100 150 200 -1 -0.5 0 0.5 1 (b) Spannungsrelaxation (blau) auf sprunghafte Deformation (rot). Abbildung 2.6: Skizzierung des modellierbaren Kriechund Relaxationsvehaltens bei Anwendung des Anti-Zener-Modells. Beliebige Zeiteinheit Spannung bzw. Verzerrung 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 (a) Kriechen (blau) durch sprunghafte mechanische Spannung (rot). Beliebige Zeiteinheit Spannung bzw. Verzerrung 0 50 100 150 200 -0.2 0 0.2 0.4 0.6 0.8 1 (b) Spannungsrelaxation (blau) auf sprunghafte Deformation (rot). Abbildung 2.7: Skizzierung des modellierbaren Kriechund Relaxationsverhaltens bei Anwendung des Zener-Modells. 2.3 Grundlagen der linearen Viskoelastizität 17 wobei τε=τσ( 1 +kλγd) gilt [Zen48]. Daraus folgt, dass sich ein ideal elastischer Körper durch Berücksichtigung der thermoelastischen Kopplung und damit innerer Reibung (Wärmeentwicklung) in einen anelastischen Körper im Sinne des Zener-Modells wandelt. Eine verallgemeinerte Betrachtung der Wärmeleitungsgleichung [Mai10; Pov11] und der adiabatischen Dehnung [Mai10] mit Hilfe fraktionaler Ableitungen führt zu einer Generalisierung der Relaxationsprozesse und damit auch zu einer Verallgemeinerung des Zener-Modells, welches als fraktionales Zener-Modell bekannt ist. Die fraktionalen Ableitungen mit der fraktionalen Ableitungsordnung ν∈R+ seien im Folgenden nach Caputo definiert [Cap67], siehe Anhang A.1. Mit der verallgemeinerten Darstellung der Temperaturänderungsmechanismen dν dtν∆T=−1 τν∆T−γdν dtνεund 0 < ν ≤1 (2.38) ergibt sich die Darstellung des fraktionalen Zener-Modells zu σ+τν σ dν dtνσ=E ε+τν ε dν dtνε!, (2.39) wodurch eine Generalisierung der im Zener-Modell modellierten exponentiellen Relaxationsund Retardationsprozesse erreicht werden kann. Anstelle der Exponentialfunktionen im Relaxationsmodul und der Kriechnachgiebigkeit, vergleiche Gl. (2.34), werden nun die Relaxationsund Retardationsprozesse durch die Mittag-Leffler-Funktion Eν(−tν) beschrieben. Die Mittag-LefflerFunktion beschreibt dabei einen Übergang zwischen hyperpolischen Funktionen ( ν= 0) und der Exponentialfunktion ( ν= 1), siehe Anhang A.2. Die Materialfunktionen des fraktionalen Zener-Modells schreiben sich nun [Mai10] J(t) = τν σ Eτν ε +τν ε−τν σ Eτν ε 1−Eν "−t τε#ν!! (2.40) G(t) = E+Eτν ε−τν σ τν σ Eν "−t τσ#ν!. Es konnte bereits gezeigt werden, dass das fraktionale Zener-Modell passend zur Beschreibung einer großen Breite von Materialien ist [Pri96; Pri03; HN11; NH13]. Die allgemeine in der Literatur beschriebene Darstellungsweise des fraktionalen Zener-Modells ist gegeben durch [Mai10; BT86; HN11; NH13] σ+τν1 σ ∂ν1 ∂tν1σ=E ε+τν2 ε ∂ν2 ∂tν2ε!, (2.41) mit ν1 , ν2∈ { 0..1 } ⊂ R+ . Bageley et al. konnten jedoch zeigen, dass das Materialmodell sich nur dann thermodynamisch sinnvoll verhält (keine negative Energiedissipationsrate, keine negative interne Arbeit), wenn folgende Bedingungen für die Parameter des Modells erfüllt sind [BT86]: E≥0 ; τν2 ε≥τν1 σ>0 ; ν1=ν2=ν(2.42) 18 2 Grundlagen viskoelastischer Materialien Somit lässt sich Gl. (2.41) wieder in Gl. (2.39) überführen. Diese Form soll im Folgenden weiter verwendet werden. Das fraktionale Zener-Modell kann auch über eine Verteilung vieler exponentieller Relaxationsvorgänge interpretiert werden [NH13; Cha05]. Es wird angenommen, dass auf mikroskopischer Ebene eine Vielzahl verschiedener Relaxationsvorgänge mit exponentiellem Charakter ablaufen. Makroskopisch kann jedoch nur die Überlagerung betrachtet werden [Cha05], welche sich dann in Superposition der Mittag-Leffler-Funktion annähern lässt. Der Relaxationsmodul des sogenannten generalisierten Zener-Modells, welcher durch eine Parallelschaltung von nZM Zener-Modellen charakterisiert ist (auch als Prony-Reihe bekannt), ergibt sich dann zu G(t) = E+ nZM X n Gne−t τσ,nmit Gn=Eτε,n−τσ,n τσ,n . (2.43) Spektrale Betrachtung ausgewählter Materialmodelle Durch Fourier-Transformation der Materialdifferentialgleichungen und Bildung des Verhältnisses zwischen komplexen mechanischen Spannungen und Deformationen lässt sich direkt der komplexe Modul G+( i ω) berechnen. Entsprechend lässt sich das modellierbare spektrale Verhalten in Anlehnung an die Messungen aus Abb. 2.1 betrachten. Das Anti-Zener-Modell schreibt sich im Fourier-Bereich GAZ +(iω) = ˜σ(iω) ˜ε(iω)=˜ηiω−ω2τε 1+iωτσ . (2.44) Abb. 2.8 zeigt exemplarisch den Verlauf des Speichermoduls und des Verlustfaktors über der Frequenz. Auffällig ist, dass für kleine Frequenzen der Speichermodul gegen Null strebt und für hohe Frequenzen einen stationären Endwert erreicht. Diese Beobachtungen decken sich mit den im Zeitbereich getroffenen Aussagen zum Relaxationsmodul G(t) . Der Verlustfaktor ist für niederfrequente Belastungen hoch, sinkt dann auf ein Minimum bei ω√τετσ= 1 und steigt dann für hohe Frequenzen wieder an. Die Frequenz des kleinsten Verlustfaktors beschreibt ebenfalls den Wendepunkt im Verlauf des Speichermoduls. In Abschnitt 4.2.1 wird das Anti-Zener-Modell noch einmal re-interpretiert, sodass an dieser Stelle diese Betrachtung zunächst ausreichen soll. Die Betrachtung der fraktionalen Ableitung im Spektralbereich vereinfacht die Anschauung erheblich. Es gilt [Cap67] dν dtνx(t)(iω)ν˜ X(iω). (2.45) Damit ergibt sich für das fraktionale Zener-Modell der komplexe Modul zu GFZ +(iω) = ˜σ(iω) ˜ε(iω)=E1+(iωτε)ν 1+(iωτσ)ν(2.46) 2.3 Grundlagen der linearen Viskoelastizität 19 Frequenz f/Hz Verlustfaktor tan δ/rad 100102104106 1 1.2 1.4 1.6 (a) Verlustfaktor. Frequenz f/Hz Speichermodul G0/GPa 100102104106 0 0.5 1 1.5 2 (b) Speichermodul. Abbildung 2.8: Beispiel der spektralen Darstellung des Anti-Zener-Modells. Parametrisierung: ˜η= 10 kPa s , τε=10 µs, τσ=40 µs. Abb. 2.9 zeigt exemplarisch für drei verschiedene fraktionale Ableitungsordnungen ν den Verlauf des Speichermoduls und des Verlustfaktors über der Frequenz. Entgegen dem Verhalten des AntiZener-Modells wird für kleine Frequenzen ein Grenzwert im spektralen Verlauf des Speichermoduls erreicht, der mit dem quasistatischen Elastizitätsmodul (wie es bspw. über Zugversuche ermittelt wird) korrespondiert und auch als relaxierter Modul bezeichnet wird. In der Parametrisierung entspricht dieser Wert dem Parameter E . Für hohe Frequenzen wird wiederum ein Grenzwert erreicht, welcher sich zu lim ω→∞Re nGFZ +(iω)o=G0 FZ,∞=E(τε/τσ)ν(2.47) berechnen lässt und als unrelaxiertes Modul bezeichnet wird. Das Maximum des Verlustfaktors liegt an der Stelle ωδ√τετσ≈ 1 und ist unabhängig von ν . Der Wendepunkt im Verlauf des Speichermoduls fällt für den Fall ν= 1 mit der Position des maximalen Verlustfaktors zusammen. Bei Variation von ν wird der Wendepunkt deutlich verschoben, das Maximum des Verlustfaktors hingegen nicht. Es ist des Weiteren auffällig, dass bezüglich des Anti-Zener-Modells der Verlauf des Verlustfaktors nun invertiert ist. Eine intuitivere Darstellung wird erreicht, indem τε durch G0 FZ,∞ ausgedrückt wird. Gl. (2.46) schreibt sich dann GFZ +(iω) = E+G0 FZ,∞(iωτσ)ν 1+(iωτσ)ν=E1+dE(iωτσ)ν 1+(iωτσ)ν. (2.48) Wird nun die fraktionale Ableitungsordnung ν variiert und das Verhältnis dE=G0 FZ,∞E−1 konstant gehalten, ergibt sich Abb. 2.10. Es ist zu erkennen, dass die fraktionale Ableitungsordnung die Steigung des Speichermoduls deutlich beeinflusst. Bei Nutzung des Zener-Modells kann die Steigung 20 2 Grundlagen viskoelastischer Materialien nicht beeinflusst werden, was häufig zu Problemen bei der Identifikation von Materialien über sehr große Frequenzbereiche führt [Pri96]. Die Lage des Maximums des Verlustfaktors verschiebt sich entsprechend der obigen Abschätzung mit ωδτσd1/(2ν) E≈ 1. Des Weiteren führt die Vergrößerung der Spannungsrelaxationskonstante τσ zu einer Verschiebung der Verläufe des Speichermoduls und des Verlustfaktors hin zu niedrigeren Frequenzen, vergleiche Abb. 2.11. In den Abbildungen nicht zu erkennen, jedoch aus den Gleichungen zu erschließen, ist dass eine Erhöhung des Verhältnisses dE zu einer Vergrößerung des Maximalwerts des Verlustfaktors führt ohne jedoch dessen Lage bezüglich der Frequenz zu verändern. Validierung der Kausalität nach Kramers-Kronig Die Kramers-Kronig Beziehungen (auch Plemelj-Gleichungen) beschreiben die Beziehung zwischen Realund Imaginärteil eines Frequenzgangs eines beliebigen Systems unter der Voraussetzung, dass das System linear und kausal ist sowie nur reell-wertige Zeitsignale zugelassen sind. Angewendet auf den komplexen Modul G+( i ω) lassen sich die Kramers-Kronig Beziehungen dann schreiben zu [Nus72; WMM05; Geo07] G0(ω) = 1 π− +∞ w −∞ G00(ω0) ω0−ωdω0(2.49) G00(ω) = −1 π− +∞ w −∞ G0(ω0) ω0−ωdω0. (2.50) In dieser Darstellung ist die Interpretation einfach über die Hilbert-Transformation gegeben. Der Realteil des Frequenzgangs ist die Hilbert-Transformierte des Imaginärteils und vice versa. Analog gelten diese Gleichungen auch für die komplexe Nachgiebigkeit J+( i ω) , allerdings mit gedrehtem Vorzeichen. Der Tausch der Vorzeichen lässt sich über die kompakte Darstellung der obigen Gleichungen veranschaulichen [Nus72]: G+(iω) = 1 iπ− +∞ w −∞ G+(iω0) ω0−ωdω0. (2.51) Gl. (2.51) und damit implizit Gl. (2.49) und Gl. (2.50) lassen sich in dieser Form nur anwenden, wenn G+( i ω) quadratisch integrierbar ist. Ist diese Bedingung nicht erfüllt, kann Gl. (2.51) unter der schwächeren Bedingung, dass der komplexe Modul eine beschränkte Funktion darstellt mit |G+(iω)|2≤C∈R+<∞, abgewandelt werden zu [Nus72] G+(iω)−G+(∞) = 1 iπ− +∞ w −∞ G+(iω0)−G+(∞) ω0−ωdω0. (2.52) 2.3 Grundlagen der linearen Viskoelastizität 21 Frequenz f/Hz Verlustfaktor tan δ/rad 100102104106 0 0.1 0.2 0.3 0.4 0.5 0.6 (a) Verlustfaktor bei Variation der fraktionalen Ableitungsordnung. Frequenz f/Hz Speichermodul G0/GPa 100102104106 0 1 2 3 4 5 (b) Speichermodul bei Variation der fraktionalen Ableitungsordnung. Abbildung 2.9: Spektrale Darstellung des fraktionalen Zener-Modells mit ν= 1 (blau), ν= 0.8 (grün) und ν=0.6 (rot). Parametrisierung: E=1 GPa, τε=40 µs, τσ=10 µs. 100102104106108 0 0.2 0.4 0.6 0.8 1 1.2 Frequenz f/Hz Verlustfaktor tan δ/rad (a) Verlustfaktor bei Variation der fraktionalen Ableitungsordnung. 100102104106 0 10 20 30 40 50 Frequenz f/Hz Speichermodul G0/GPa (b) Speichermodul bei Variation der fraktionalen Ableitungsordnung. Abbildung 2.10: Spektrale Darstellung des fraktionalen Zener-Modells mit ν= 1 (blau), ν= 0.8 (grün) und ν=0.6 (rot). Parametrisierung: E=1 GPa, G0 FZ,∞=50 GPa, τσ=10 µs. 100102104106108 0 0.2 0.4 0.6 0.8 1 1.2 Frequenz f/Hz Verlustfaktor tan δ/rad (a) Verlustfaktor bei Variation der Spannungsrelaxationskonstante. 100102104106 0 10 20 30 40 50 Frequenz f/Hz Speichermodul G0/GPa (b) Speichermodul bei Variation der Spannungsrelaxationskonstante. Abbildung 2.11: Spektrale Darstellung des fraktionalen Zener-Modells mit τσ= 1 µs (blau), τσ= 10 µs (grün) und τσ=100 µs (rot). Parametrisierung: E=1 GPa, G0 FZ,∞=50 GPa, ν=1. 28 3 Stand der Forschung zur Materialparameterbestimmung und akustischer Wellenausbreitung Einträge in der Steifigkeitsmatrix zugelassen. Zuverlässige Ergebnisse wurden bezüglich der Realteile der E-Moduln erzielt, jedoch wurde eine erhöhte Unsicherheit in Bezug auf deren Imaginärteile beobachtet. Ein konkretes viskoelastisches Modell wurde nicht verwendet, sodass die Kausalität dieses Modells in Frage zu stellen ist. Ein ähnlicher Ansatz wurde Anfang der 1990er Jahre von Sachse et al. unter dem Namen Point-Source/Point-Receiver Methode vorgestellt [ES90; CKS91]. Hier wird mit einem Laser räumlich stark begrenzt eine mechanische Welle erzeugt und mit einem Empfänger schmaler Apertur auf der gegenüberliegenden Seite des Prüfobjekts unter verschiedenen Winkeln bezogen auf den direkten Weg vom Sender zum Empfänger abgetastet. Durch die impulsförmige Anregung werden eine Quasi-Longitudinalund zwei Quasi-Transversalwellen angeregt, deren Eintreffzeitpunkte am Empfänger separiert werden müssen. Mit Hilfe eines inversen Ansatzes kann das Materialmodell (reell-wertige Moduln) identifiziert werden. Als Modell wurde wiederum die Christoffel-Gleichung verwendet und die quadratische Abweichung der berechneten und experimentell bestimmten polaren Winkel der ebenen Wellenfronten durch einen Quasi-Newton Algorithmus minimiert. Probleme der Methode bestehen in der Identifikation bzw. Separation der verschiedenen Wellenarten und der Diskrepanz zwischen gemessener Gruppenund für das Modell relevanter Phasengeschwindigkeit. Eine Laufzeitmessung spiegelt immer eine Gruppenlaufzeit wider und kann streng genommen nur unter Freiraumbedingungen und Vernachlässigung von Absorptionsphänomenen als äquivalent zur Phasengeschwindigkeit angesehen werden. Letzteres Problem gilt auch für die Immersionstechnik, wobei bei der Immersionstechnik zusätzlich problematisiert wird, dass sich viele Proben durch das Eintauchen in das Wasserbad hinsichtlich ihrer Eigenschaften verändern. Eine bezüglich der Umgebungseinflüsse konditionierte Messung ist somit ausgeschlossen. Abhilfe schafft die Luft-gekoppelte Anregung von Wellen [SLC96; SFN10], welche jedoch aufgrund schlechter Impedanzanpassung hinsichtlich des erzielbaren Signal-Stör-Abstands (SNR) und des eingeschränkten Frequenzbereichs aufgrund der hohen Güte von Luftschallwandlern problematisch ist. Geführte Wellen Ansätze Parallel zu den Ansätzen mit ebenen Wellen, wurden in den letzten Jahren auch geführte Wellen Ansätze zur Bestimmung der Materialparameter von Probekörpern entwickelt. Anders als bei den ebene Wellen Ansätzen wird hier gezielt die räumlich beschränkte Ausdehnung der Proben konstruktiv genutzt. Entsprechend ist die geometrische Form der Proben für das hinterlegte Modell bedeutsam. Meist werden plattenförmige Probekörper genutzt, um das dispersive Verhalten der symmetrischen und antisymmetrischen Lamb-Wellen auszuwerten. Im Wesentlichen unterscheiden sich geführte Wellen Ansätze durch die gewählte Messmethode, d.h. den instrumentellen Einsatz, zur Bestimmung der dispersiven Kenngrößen des Wellenleiters, durch die gewählte Modellierungsmethode der geführten Wellenausbreitung und durch das zu identifizierende Materialmodell. 3.1 Materialparameterbestimmung 29 Vishnuvardhan et al. zeigten im Rahmen einer Sensitivitätsstudie, dass die Ausbreitungsgeschwindigkeit der symmetrischen Grundmode in einem orthotropen plattenförmigen Material kaum sensitiv ist, die antisymmetrische Grundmode dagegen deutlich auf die E-Moduln C44 und C55 reagiert [VKB07b]. Es ist demnach sinnvoll, eine multimodale Auswertung durchzuführen. Identifiziert wurde ein orthotroper elastischer Modellansatz mit Hilfe eines inversen Ansatzes ausgeführt über einen genetischen Optimierungsalgorithmus, wobei ein analytisches Modell (Rayleigh-Lamb Gleichung) genutzt wurde. Von messtechnischer Seite wurde eine Piezokeramik als Sender zentral auf eine Platte geklebt und mehrere Empfänger um diese herum angebracht. Als Empfänger dienten ebenfalls Piezokeramiken und zum Vergleich ein Laser-Doppler-Vibrometer. Instrumentativ deutlich weniger aufwendig nutzt Rogers beispielsweise einen Wedge-Transducer, welchen er entlang eines plattenförmigen Wellenleiters verschiebt, um über einen Phasenvergleich mit dem Sendesignal die Randwellenlänge der Leaky-Lamb-Welle zu bestimmen [Rog95]. Rogers bestimmt mit diesem Messaufbau ein isotropes elastisches Materialmodell. Instrumentativ aufwendiger ist die räumlich äquidistante Abtastung des zeitlichen Verlaufs der Wellenbewegung am Wellenleiterrand. Durch Berechnung einer 2D-Fourier-Transformation der räumlich und zeitlich aufgelösten Daten wird das (experimentelle) Dispersionsdiagramm der angeregten (und am Wellenleiterrand sichtbaren) Wellenleitermoden 3 erstellt. Dabei wird aus der Zeitachse die Frequenzachse und aus der Raumachse die Achse der Wellenzahlen. Die Reproduzierbarkeit der Messung und der Anregung ist für diese Methode von Bedeutung, sodass häufig kontaktlose Anregungen und Messungen über eine hoch-energetische Laserquelle und ein Laser-Doppler-Vibrometer realisiert werden, vgl. z.B. [HJQ00; Sin15]. Zur Reduktion des instrumentellen Aufwands kann z.B. das kostenintensive Laser-Doppler-Vibrometer durch einen Miniatur-Transducer ersetzt werden [HC12], welcher mittels einer Lineareinheit verschoben wird. Eine Forschergruppe um Marzani nutzt hingegen zur weiteren Reduzierung des technischen Aufwands und zur Umgehung des Problems der schwachen Reproduzierbarkeit der Einzelmessungen nur einen örtlich festen piezokeramischen Empfangswandler und eine örtlich feste Laserquelle und berechnet aus dem aufgezeichneten Zeitsignal auf Basis der Gabor-Wavelet-Transformation den frequenzabhängigen Verlauf der Gruppengeschwindigkeit der symmetrischen und antisymmetrischen Grundmode [SRM11]. Die Arbeitsgruppe geht noch einen Schritt weiter und substituiert die Laserquelle durch einen piezoelektrischen Schallwandler [MM12]. In [SRM11] wird ein isotropes elastisches Materialmodell mit Hilfe des Nelder-Mead Simplex-Algorithmus invers charakterisiert. Hingegen wird in [MM12] ein orthotropes elastisches Materialmodell identifiziert und für den inversen Ansatz ein genetischer Optimierungsalgorithmus angesetzt. Identisch ist hingegen der Ansatz, aus den Zeitsignalen über Zeit-Frequenz-Analyse die frequenzabhängigen Gruppengeschwindigkeiten der angeregten Wellenleitermoden zu berechnen und mit simulierten Daten zu vergleichen. Als Simulationsmodell wird in beiden Publikationen die SAFE Methode, vgl. nächster Abschnitt, verwendet. Ein erweiterter Überblick über geführte 3 Im Folgenden werden Wellenleitermoden auch kurz Moden genannt, solange im Kontext die Eindeutigkeit sichergestellt ist. 30 3 Stand der Forschung zur Materialparameterbestimmung und akustischer Wellenausbreitung Wellen in Plattenstrukturen sowie deren Möglichkeiten in der Charakterisierung von Materialien ist in [Chi97] gegeben. Im Gegensatz zu plattenförmigen Wellenleiterstrukturen ist die Literaturlage bezüglich der Materialcharakterisierung mittels zylindrischer Probekörper deutlich dünner. Dies liegt wohl nicht zuletzt an den zahlreichen Arbeiten bezüglich der Charakterisierung von gewalzten Metallplatten und zeitlich später den Komposit-Platten. Während Rohrgeometrien bei der Charakterisierung von (viskoelastischen) Fluiden eine zentrale Rolle in der ultraschallbasierten Prozessmesstechnik spielen [SC03; VLC04] 4 , wird die zylindrische Struktur zur Materialcharakterisierung von Polymeren bislang wenig verfolgt. Extrudierte Stäbe (und auch Platten) spielen wirtschaftlich eine große Rolle, wenn aus technischen oder wirtschaftlichen Gründen kein Spritzgussverfahren eingesetzt werden kann. Die extrudierten Halbzeuge können nach dem Abkühlen spanend bearbeitet oder auch warm umgeformt werden. Gerade im Falle kleiner Stückzahlen, komplexer Geometrien oder großer Wandstärken wird dieses Verfahren häufig favorisiert. Einen ersten Ansatz zur Messung von Longitudinalund Transversalwellengeschwindigkeit in einem zylindrischen Probekörper liefert Reynolds im Jahr 1953 [Rey53]. Der Durchmesser des Stabs wird als groß gegenüber der betrachteten Wellenlänge beschrieben, sodass ein Strahlenansatz ohne Berücksichtigung von Wellenleitermoden angesetzt werden kann. Identifiziert wird in einem direkten Verfahren ein isotropes elastisches Modell von metallischen Proben. Rautenberg griff diesen Ansatz 2012 wieder auf, um stark absorbierende Kunststoffe zu charakterisieren [Rau12]. Die Probekörper wurden deutlich verkleinert, sodass wellentheoretische Ansätze zur Lösung des Wellenleiterproblems notwendig wurden, jedoch konnte der Ansatz von Reynolds zur Gewinnung von Startwerten für das inverse Problem, welches mit dem Nelder-Mead Simplex-Algorithmus gelöst wurde, verwendet werden. Identifiziert wurde von Rautenberg ein elastisches isotropes und ein transversal isotropes (mit quasi-isotroper Näherung) Materialmodell sowie zwei Dämpfungsparameter, welche zum Rayleigh-Dämpfungsmodell gehören. Die Arbeiten von Rautenberg bilden die Grundlage für die vorliegende Arbeit. Übernommen werden in zentralen Punkten die Probengeometrie, das Konzept der Startwertbestimmung mit einem Strahlenansatz nach Reynolds [Rey53] sowie das Konzept des inversen Verfahrens mit einem wellentheoretischen Modell zur Simulation geführter Wellen. Zentrale Abgrenzungspunkte zur vorangegangenen Arbeit von Rautenberg sind •das zu identifizierende Materialmodell, •der Lösungsansatz des Vorwärtsmodells, •der betrachtete Frequenzbereich, •der Aufbau, die Modellierung und die Identifikation der Schallwandler sowie •die Betrachtung von Unsicherheiten des Messsystems. 4Vergleiche z.B. auch die Liner Serie der SensAction AG. URL: www.sensaction.de 3.2 Modellierung geführter Wellenausbreitung 31 3.2 Modellierung geführter Wellenausbreitung Im Laufe der letzten Jahrzehnte wurden viele Modellierungsansätze zur Analyse von Wellenausbreitungsphänomenen in Wellenleiterstrukturen entwickelt und erfolgreich angewendet. Aktuell konkurrieren im Wesentlichen vier unterschiedliche Methoden in der wissenschaftlichen Szene: Die Global Matrix Method (GMM) [Low95], die Waveguide FEM (WFEM) [MDB05], die SemiAnalytic Finite Element Method (SAFE) [MVB08] und die Scaled Boundary Finite Element Method (SBFEM) [Gra14; GBS14]. Die GMM verfolgt einen analytischen Ansatz mittels analytischer Formulierung von geometrieund materialbedingten Ansatzfunktionen. Es wird ein lineares homogenes Gleichungssystem aufgestellt, dessen nicht-triviale Lösungen die ausbreitungsfähigen Wellenleitermoden charakterisieren. Problematisch ist in diesem Zusammenhang die numerische Nullstellensuche in der Determinanten der Systemmatrix. Häufig werden nicht alle Nullstellen der Funktion gefunden und somit modale Lösungen ausgelassen. Ein Ansatz dieses Problem zu lösen ist die Berechnung aller Nullstellen mit Hilfe des Intervall-Newton-Verfahrens [WBH12a; BWH12; BWR13; BUW13], welches sich jedoch, aufgrund langer Rechenzeit, nicht zielführend innerhalb eines inversen Problems zeigt. Des Weiteren ist die Berechnung der Determinanten in der Nähe einer Nullstelle numerisch gesehen eine Herausforderung [BWH12; WBH12b]. Weiterer Kritikpunkt an der Methode ist der Übergang zu komplexen Nullstellenproblemen bei der Verwendung komplexer Materialformulierungen. Die WFEM und die SAFE Methode basieren beide auf Standard FEM Werkzeugen. Während bei der WFE Methode ein Einheitssegment des Wellenleiterquerschnitts mittels FEM diskretisiert wird, benötigt die SAFE Methode nur den tatsächlichen Querschnitt des Wellenleiters. Die SAFE Methode kann für plattenförmige oder zylindrische Wellenleiter auch auf eine diskretisierte Dimension reduziert werden [MVB08]. In beiden Methoden werden die FE-Matrizen restrukturiert, sodass sich ein Eigenwertproblem zur Lösung des jeweiligen Systems formulieren lässt. Werden komplexe Materialmodelle angenommen, werden die Eigenwerte komplex, wobei die Lösung dieses Problems numerisch gesehen deutlich einfacher zu handhaben ist als eine Nullstellensuche in der komplexen Ebene. Durch die Diskretisierung eines Wellenleiterquerschnittsegments bei der WFE Methode wird die auf Eigenwerte und Eigenvektoren zu analysierende Matrix sehr groß, was sich negativ auf die Rechenzeit und auch auf die Stabilität der Lösung auswirkt. Da für Platten und Zylinder eine Reduktion der Diskretisierung auf nur eine Dimension ausreichend ist, um die geführte modale Wellenausbreitung zu betrachten, ist die SAFE Methode bezüglich ihrer Stabilität und vor allem der Rechengeschwindigkeit deutlich zu favorisieren. Eine noch junge Methode zur Berechnung von geführten akustischen Wellen ist die SBFEM. Grundlegend wurde die SBFEM für elastodynamische Probleme von Song 2009 beschrieben [Son09]. Darauf aufbauend folgte eine Reihe wissenschaftlicher Publikationen von Gravenkamp, welche sich speziell mit der Ausbreitung geführter elastischer Wellen im Kontext der SBFEM beschäftigen. So entstanden Arbeiten zur Modellierung von Plattenwellenleitern [GSP12; Gra14], zylindrischen 32 3 Stand der Forschung zur Materialparameterbestimmung und akustischer Wellenausbreitung Wellenleitern [GBS14a; GBS14; BGS14; Gra14], Wellenleitern beliebigen Querschnitts [GMS13; Gra14] sowie eingebetteten Wellenleitern [GBS14b; Gra14]. Grundlegende Idee der SBFEM ist die Diskretisierung eines Rechengebietes mittels Finiter Elemente ausschließlich an dessen Rand. Beispielsweise wird ein zweidimensionales Rechengebiet mit Linienelementen am Rand diskretisiert und ein Skalierungszentrum eingeführt, mit welchem der Rand kontinuierlich skaliert werden kann und somit jeder Punkt im Raum beschrieben ist. Dies führt dazu, dass sich im Vergleich zur FEM die Dimensionalität der Formfunktionen der Finiten Elemente um eine Ordnung reduziert und die Knotenpunktverschiebungen gemischt (an den Knotenpunkten der Finiten Elemente) diskret und analytisch (kontinuierlich) mit der Skalierungsvariablen vorliegen. Wesentliche Vorteile dieser Methode hinsichtlich der Modellierung geführter akustischer Wellen gegenüber der SAFE Methode sind niedrigere Rechenzeiten (gerade bei hohen Frequenzen) und eine einfachere Möglichkeit zur Beschreibung angrenzender Medien an den Wellenleiter [Gra14]. Der Rechenvorteil gegenüber der SAFE Methode entsteht durch kleinere Matrizen zur Berechnung des (einfachen) Eigenwertproblems sowie der Hamiltonischen Eigenschaft der Matrizen, welche spezielle Lösungsalgorithmen mit hoher Effizienz zulassen. 3.3 Transiente Finite Elemente Simulation in viskoelastischen Medien Dieser Abschnitt charakterisiert den aktuellen Stand der Zeitbereichssimulation in der FEM hinsichtlich der Implementierbarkeit viskoelastischer Verlustmechanismen. Für die vorliegende Arbeit sind die aus dieser Entwicklung ableitbaren Konsequenzen für die Materialparameterbestimmung von Bedeutung. Der Stand der Technik wird im Folgenden grob umrissen. Es wird im Kontext dieser Arbeit kein neues numerisches Verfahren zur Zeitbereichssimulation in der FEM vorgelegt. Die Finite Elemente Methode (FEM) ist ein leistungsfähiges numerisches Lösungsverfahren partieller, zeitund ortsabhängiger Differentialgleichungen mit definierten Anfangsund Randbedingungen und hat sich zu einem Standard-Verfahren für die Berechnung physikalischer Systeme entwickelt [Kal04]. Vorteile der FEM liegen aus numerischer Sicht in ihrer Effizienz. Die Systemmatrizen sind zumeist symmetrisch und dünn besetzt, was die Speicherung und Lösung der Gleichungssysteme begünstigt. Aus Anwendersicht liegt der Vorteil in der Möglichkeit auch komplexe Geometrien zu analysieren, wobei sowohl statische, transiente, harmonische Analysen als auch Eigenfrequenzanalysen möglich sind [Kal04]. Zwei wesentliche Grenzen der Methode seien gerade hinsichtlich der Simulation von Ultraschallwellen genannt. Zum einen ist der räumliche Diskretisierungsaufwand des Rechengebiets sehr schnell groß und zum anderen sind ’offene’ Randbedingungen zur Simulation eines unendlich ausgedehnten Raums problematisch [Kal04]. 3.3 Transiente Finite Elemente Simulation in viskoelastischen Medien 33 3.3.1 Grundgleichungen und das Rayleigh-Dämpfungsmodell in der Finite Elemente Methode Ausgangspunkt der Betrachtung für elastodynamische Probleme ist die Cauchysche Bewegungsgleichung. In der hier dargelegten Form sind bereits die Annahmen eines homogenen Hooke’schen Materials, d.h. linear elastisches Materialverhalten, enthalten. Mit dem Volumenkraftvektor fV und der Steifigkeitsmatrix Cergibt sich [Kal04] ΥTCΥξ+fV=ρ∂2ξ ∂t2, (3.1) wobei Υeine Matrix partieller Ableitungen beschreibt ΥT= ∂ ∂x0 0 0 ∂ ∂z ∂ ∂y 0∂ ∂y0∂ ∂z0∂ ∂x 0 0 ∂ ∂z ∂ ∂y ∂ ∂x0  . (3.2) Zunächst wird die Cauchysche Bewegungsgleichung nach Gleichung (3.1) in ihre schwache Form überführt. Dazu wird eine Testfunktion ξ0 mit der partiellen DGL multipliziert und über das Rechengebiet Ω partiell integriert. Unter Verwendung des Green’schen Satzes und dem Ausblenden der Beiträge vom Rand ∂Ω über entsprechende Wahl der Randbedingungen (Details siehe [Kal04]) ergibt sich w Ω ρξ0∂2ξ ∂t2dΩ+w Ω (Υξ0)TCΥξdΩ−w Ω ξ0fVdΩ=0. (3.3) Eine räumliche Diskretisierung in Subräume Ωe mit Ω=PΩe soll nun den räumlich kontinuierlichen mechanischen Verschiebungsvektor ξhinreichend approximieren. Es gilt ξ≈ ne X b=1 UbNb;Nb= Nb0 0 0Nb0 0 0 Nb  , (3.4) wobei ne die Anzahl der Knotenpunkte ohne Dirichlet Randbedingungen, Ub die diskreten Knotenpunktverschiebungen am Knoten b und Nb eine passende Interpolationsfunktion (auch Formfunktion) beschreiben. Die gleiche Approximation wird auch für die Testfunktion durchgeführt. Einsetzen in die schwache Form der Cauchyschen Bewegungsgleichung liefert die semidiskrete Galerkin Formulierung ne X a=1 ne X b=1w Ω ρNT aNbdΩ·∂2Ub ∂t2+w ΩNaΥTTCNaΥTdΩ·Ub−w Ω NT afVdΩ=0, (3.5) 34 3 Stand der Forschung zur Materialparameterbestimmung und akustischer Wellenausbreitung welche in Matrixform in die bekannte Form der FE-Gleichungen umgeschrieben werden kann zu Mu¨ U+KuU=F. (3.6) Die Matrizen werden als Massenmatrix Mu und Steifigkeitsmatrix Ku bezeichnet. U bezeichnet die diskreten unbekannten Verschiebungen an den Knotenpunkten. Das Gleichungssystem beschreibt formal ein Feder-Masse-System mit Kraftanregung, d.h. ein ungedämpftes schwingungsfähiges System. Schwingungsfähige Systeme zeigen in der Realität immer ein Dämpfungsverhalten. In (konventionellen; kommerziell verfügbaren) FEM-Werkzeugen werden diese Verluste durch eine zusätzliche Dämpfungsmatrix Cu berücksichtigt, welche in Anlehnung an das Kelvin-Voigt-Modell als geschwindigkeitsproportionaler Faktor angesetzt wird. Gleichung (3.6) wird somit erweitert zu Mu¨ U+Cu˙ U+KuU=F, (3.7) wobei sich Cu als Linearkombination aus Massenund Steifigkeitsmatrix mit den Koeffizienten αM und αKdarstellen lässt: Cu=αMMu+αKKu. (3.8) Damit bleiben formal alle Matrizen reell-wertig, obwohl nun Dämpfung eingebunden worden ist. Wichtig ist zu erwähnen, dass die Struktur des Rayleigh-Dämpfungsmodells keiner physikalischen Interpretation entspringt, sondern ein mathematisches Konstrukt ist, welches äquivalent wirkende Dämpfungseffekte mit möglichst wenigen zusätzlichen Parametern beschreibt. Außerdem lässt sich diese Art der Dämpfungsformulierung leicht in Zeitschrittverfahren einbinden. Die ZeitschrittDiskretisierung hyperbolischer partieller DGL wird häufig mit dem Newmark-Verfahren realisiert [Kal04]. Es gehört zu den sogenannten Einschritt-Verfahren, was bedeutet, dass nur ein vergangener zeitlicher Zustand für die Berechnung des nächsten Zeitschritts gespeichert werden muss. Wie von Rautenberg [Rau12] gezeigt, treten bei Verwendung des Rayleigh-Dämpfungsmodells nicht-kausale Effekte auf. Rautenberg untersuchte dazu u.a. empirisch den Einfluss der RayleighDämpfung. So modellierte er eine akustische Übertragungsstrecke definierter Länge, welche im Zeitschrittverfahren durch transiente FE-Simulationen bei verschiedenen Dämpfungsparametern simuliert wurde. Es konnte beobachtet werden, dass sich die simulierten Eintreffzeitpunkte des jeweiligen Wellenpakets in Abhängigkeit von der eingestellten Schallabsorption verschieben. Eine höhere Schallabsorption des Materials führte dabei zu einer Verschiebung hin zu früheren Eintreffzeitpunkten. Grundlegend ist das Problem der Rayleigh-Dämpfung darin zu sehen, dass im äquivalenten Sinne ein frequenzabhängiges Dämpfungselement im Kelvin-Voigt-Modell angesetzt wird (vergleiche Abschnitt 4.2.1). Dieser Ansatz ist in der Literatur bezüglich Kausalitätsproblemen schon früher und losgelöst von der konkreten Implementierung in der FEM diskutiert worden [Cra70; Adh13]. 3.3 Transiente Finite Elemente Simulation in viskoelastischen Medien 35 3.3.2 Zeit-kausale viskoelastische Dämpfung in der Finite Elemente Methode Wie in Kapitel 2.3 erläutert, verändert sich die Beziehung zwischen mechanischen Spannungen und Dehnungen durch die Berücksichtigung des zeitlichen Beanspruchungsgedächtnisses viskoelastischer Materialien. Nach Gl. (2.15) ergibt sich, nun für den dreidimensionalen Fall notiert, ein Faltungsintegral zwischen der zeitlichen Ableitung des mechanischen Verzerrungstensors [˙] und einem zeitabhängigen Tensor 4-ter Stufe [G(t)] , welcher die Relaxation des Werkstoffs beschreibt. σi j(t) = Gi jkl(t)∗˙kl(t)(3.9) Wird diese Gleichung in die Bewegungsgleichung von Cauchy (3.1) eingesetzt, sieht man leicht, dass sich die Faltung bis in die semidiskrete Galerkin Formulierung fortsetzt. Formal ergibt sich ein Ausdruck der Form [SL01] ˜ Mu·¨ U+˜ Ku(t)∗˙ U=˜ F. (3.10) Limitierend für die FE-Simulation ist an dieser Stelle, dass ohne Umformulierung dieses Problems zu jedem Knotenpunkt die volle Historie der Verschiebungen gespeichert werden muss, wenn eine Lösung im Zeitintegrationsverfahren angestrebt wird. Gerade hinsichtlich des hohen Diskretisierungsaufwands bei der Ultraschallsimulation führt dies schnell zu unhandhabbaren Datenvolumina. Aktuelle Forschungsarbeiten zeigen jedoch Ausbzw. Umwege aus diesem Dilemma. Stucky und Lord [SL01] modellieren in der FEM ein lineares isotropes Material mit dem ZenerModell, d.h. die Relaxation ist im Zeitbereich über eine Exponentialfunktion beschrieben. Es können für Scherund Volumeneigenschaften unterschiedliche Relaxationszeiten angenommen werden. Durch diesen Ansatz ist es möglich, das Faltungsintegral durch Speicherung von nur zwei vergangenen Zeitschritten auszuführen. Im Vergleich zum Newmark-Verfahren mit dem RayleighDämpfungsmodell ergibt dies jedoch trotzdem einen doppelt so großen Speicherbedarf. Des Weiteren wird auch an FE-Implementierungen für viskoelastische Wellenausbreitung basierend auf dem generalisierten Zener-Modell, vgl. Gl. (2.43), gearbeitet, wobei bereits erste Erfolge erzielt werden konnten [BEJ04; GT06; LP11]. Es ist demnach sinnvoll und konsequent, ein Messsystem zu entwickeln, welches in der Lage ist, die gängigen Modelle zur Beschreibung viskoelastischen Materialverhaltens zu identifizieren, um diese Parameter für die realitätsnahe FE-Simulation praktisch nutzbar zu machen. 4 Wellenausbreitung in viskoelastischen Medien 37 4 Wellenausbreitung in viskoelastischen Medien 4.1 Implikationen viskoelastischer Medien auf die Wellengleichungen Im Folgenden werden die Wellengleichungen in räumlich unendlich ausgedehnten linear viskoelastischen Materialien betrachtet. Dabei werden nur die im Kontext der Arbeit relevanten Grundzüge berücksichtigt. Für ausführlichere Darstellungen sei bspw. auf [CC86; CKK88] verwiesen. Ausgangspunkt der Betrachtung bildet das Materialmodell. Zur Vereinfachung der Betrachtungen sei zunächst ein homogenes und isotropes Materialverhalten angenommen. Aus den vorangegangenen Kapiteln ist bereits bekannt, dass viskoelastische Materialien ein Gedächtnis besitzen und damit die unabhängigen Komponenten der Materialbeschreibung Funktionen der Zeit darstellen 5 . Das Materialmodell -aufgrund der kurzen Schreibweise an dieser Stelle in Tensor-Notationbeschrieben über den Kompressionsmodul K(t)und den Schermodul µ(t)liefert σi j =3K(t) 1 3εkkδi j!+2µ(t) εi j −1 3εkkδi j!. (4.1) Die Bewegungsgleichung von Cauchy schreibt sich komponentenweise ρ¨ ξi=fi+ 3 X j=1 ∂σi j ∂xj . (4.2) Im Folgenden wird zur Vereinfachung der Gleichung auf die Vektordarstellung in der Voigt’schen Notation übergegangen. Durch Einbeziehung des Boltzmann’schen Superpositionsprinzips aus Gl. (2.15) in die Bewegungsgleichung unter Vernachlässigung der Startbedingung für t→ 0 + ergibt sich [CC86] ρ¨ ξ= t w 0 K(t−τ) + 1 3µ(t−τ)!∂ ∂τ∇∇·ξdτ+ t w 0 µ(t−τ)∂ ∂τ∇2ξdτ. (4.3) 5Für den isotropen Fall gültig. Diskussion bezüglich anisotroper Materialen siehe Kapitel 4.2.3. 44 4 Wellenausbreitung in viskoelastischen Medien Norm. Kreisfrequenz ω·τσ α(ω)Np/m 10−3100105 10−5 100 105 Norm. Kreisfrequenz ω·τσ cph(ω)/c0 10−3100105 1 2 3 4 (a) G0 FZ,∞= 10 GPa , ν= 1 (blau); ν= 0,8 (grün); ν=0,6 (rot). Norm. Kreisfrequenz ω·τσ α(ω)Np/m 10−3100105 10−5 100 105 Norm. Kreisfrequenz ω·τσ cph(ω)/c0 10−3100105 2 4 6 8 10 (b) G0 FZ,∞= 10 GPa (blau); G0 FZ,∞= 50 GPa (grün); G0 FZ,∞=100 GPa (rot) , ν=1. Abbildung 4.3: Darstellung des Effektes verschiedener Parametrisierungen des fraktionalen Zener-Modells auf die Phasengeschwindigkeit und die Schallabsorption. Konstante Parametrisierungen für beide Grafiken: E=1 GPa, τσ=10 µs. Es stellt sich folglich die Frage nach dem für thermoplastische Kunststoffe zu erwartenden Parameterbereich. Dieser ist im Vorfeld nur schwer abschätzbar. Einen Anhaltspunkt liefert die Einteilung der Parameter in drei Segmente und deren Betrachtung bezüglich der Frequenzabhängigkeit der akustischen Größen in diesen Bereichen. Tabelle 4.1 zeigt die Einteilung der Segmente nach [HN11]. Wird zunächst nur die Proportionalität zwischen Kreisfrequenz und Absorption betrachtet, lässt sich eine grobe Einteilung vornehmen, wobei zunächst ν= 1 gelte. Der Bereich ωτ−1 ε zeigt einen quadratischen Verlauf der Absorption mit der Frequenz, was i.A. für Schmelzen (Flüssigkeiten) und Gase angesetzt wird [Kut88]. Das überrascht hinsichtlich des Zeit-Temperatur-Superpositionsprinzips wenig, da eine geringe Frequenz einer hohen äquivalenten Temperatur entspricht und somit nach Abb. 2.1 der Schmelze-Bereich erreicht wird. Für ωτ−1 σ zeigt sich ein lineares Verhalten der Absorption bezüglich der Frequenz, was i.A. Metallen bzw. rigiden Festkörpern zugesprochen wird. Nach dem Zeit-Temperatur-Superpositionsprinzip entspricht dieser Bereich niedrigen Temperaturen im Glaszustand. Entsprechend beschreibt das mittlere Segment mit τ−1 εωτ−1 σ den Erweichungsbereich bzw. die Grenzbereiche des Glaszustands und des gummi-elastischen Bereichs, also den Bereich in der Umgebung der Glasübergangstemperatur. Aus der Beobachtung heraus, dass sich Kunststoffe bzw. Thermoplaste hinsichtlich ihrer frequenzabhängigen Absorption weder linear noch quadratisch verhalten, lässt sich im Folgenden der Anwendungsbereich auf das mittlere Segment begrenzen. 4.2 Viskoelastische Materialmodelle für die Simulation akustischer Wellen in Polymeren 45 Tabelle 4.1: Das Dispersionsverhalten des fraktionalen Zener-Modells, betrachtet in drei Segmenten nach [HN11]. Segment-Definition Phasengeschwindigkeit Absorption ωτ−1 εcph(ω)≈c0α(ω)∝ω1+ν τ−1 εωτ−1 σcph(ω)∝ων/2α(ω)∝ω1−ν/2 ωτ−1 σcph(∞) = c0τετ−1 σν/2α(ω)∝ω1−ν 4.2.3 Generalisierung auf dreidimensionale Betrachtungen und Diskussion der Anisotropie Aufgrund des Einflusses des Verarbeitungsprozesses bei der Herausbildung von charakteristischen Eigenschaften eines Polymers, werden in dieser Arbeit ausschließlich extrudierte Stäbe zur Prüfung verwendet. Ausgehend von diesem Herstellungsprozess ist es opportun anzunehmen, dass sich durch den Verarbeitungsprozess die Polymerketten vorzugsweise in Extrusionsrichtung ausrichten [CLY92; War75] und damit ein (schwach ausgeprägter) anisotroper Körper ensteht. Es stellt sich die Frage, wie nun die eindimensionalen linear viskoelastischen Modelle auf dreidimensionale anisotrope Körper übertragen werden können. Die allgemeinste Formulierung für dieses Problem kann in Anlehnung an Gl. (2.10) und Gl. (2.15) geschrieben werden zu [Chr82; Car90; DSJ89] σ(t) = ε(t=0+)Γ(t) + t w 0 Γ(t−τ)∂ε(τ) ∂τ dτ, (4.26) wobei Γ(t) eine Matrix beschreibt, deren Elemente die Relaxationsmoduln bilden. Aus dieser allgemeinen Darstellung heraus lässt sich zunächst nur ableiten, dass jedes Element der Matrix eine zeitabhängige Funktion darstellt, jedoch ist noch keine Aussage über die Zusammensetzung der zeitabhängigen Funktionen in Abhängigkeit von der Anisotropie des Körpers getroffen. Douven et al. [DSJ89] verfolgen einen direkten Weg zur Generalisierung auf die transversale Isotropie, indem sie jeden (für den elastischen Fall) linear unabhängigen elastischen Parameter im viskoelastischen Fall als zeitabhängige Funktion ansetzen. Caricone [Car90] hingegen beschreibt einen bezüglich der Anisotropie allgemeinen Ansatz, wobei er jedoch Restriktionen bezüglich der Allgemeingültigkeit des Dämpfungsverhaltens einführt. Es wird die Relaxationsmatrix Γzerlegt in drei Anteile Γ(t) = C0+C1f1(t) + C2f2(t), (4.27) wobei fi(t) mit i= 1,2 Funktionen der Zeit darstellen und in ihrer Struktur als generalisiertes Zener-Modell formuliert sind. Die Matrizen Ci werden als ideal elastisch angesetzt. Die Zerlegung der Matrix Γ wird derart gestaltet, dass Quasi-Dilatationsbewegungen des Körpers mit der Rela- 46 4 Wellenausbreitung in viskoelastischen Medien xationsfunktion f1(t) und Quasi-Scherbewegungen mit der Relaxationsfunktion f2(t) gedämpft werden. Die Zerlegung der Anteile erfolgt in der Tensor-Notation über die Entkopplung der Spur des Spannungstensors und des Spannungsdeviators. In dieser Arbeit wird ein ähnlicher, aber physikalisch verallgemeinerter Ansatz zur Generalisierung auf dreidimensionale Betrachtungen beschrieben. In Anlehnung an die spektrale Zerlegung von Materialtensoren nach Theocaris [TS01; The00; TP90] wird im Folgenden die Materialmatrix in ihre orthogonalen Anteile zerlegt. Es wird begonnen mit der elastischen Materialmatrix. Die Zerlegung ermöglicht es, im Nachhinein den orthogonalen Eigenbewegungen eines Materials viskoelastische Verluste zuzuordnen. Spektrale Zerlegung der Materialmatrix Nach Gl. (2.10) ist die elastische Steifigkeitsmatrix definiert als C∈R6×6 . Die elastische Nachgiebigkeitsmatrix S ist entsprechend als Inverse definiert mit S=C−1 . Zunächst werden alle Eigenwerte λi∈R+ und Eigenvektoren pi∈R6 der elastischen Nachgiebigkeitsmatrix berechnet. Die Eigenvektoren werden zu einer Matrix P=[p1,p2,. . .]∈R6×6 zusammengesetzt. Es gilt dann diag(λi) = P−1SP, (4.28) wobei diag(λi) eine Diagonalmatrix mit den λi auf der Hauptdiagonalen beschreibt. Das Hooke’sche Gesetz nach Gl. (2.10) kann nun umgeschrieben werden zu ε=Sσ=Pdiag(λi)P−1σ(4.29) =X i λiPdiag(δii)P−1σ(4.30) =X i λiDiσ(4.31) =X i λiσi, (4.32) wobei δii das Kronecker-Symbol beschreibt. Die orthogonalen und idempotenten Matrizen Di beschreiben eine orthogonale Projektion des Spannungsvektors σ auf die Eigen-Spannungsvektoren σi. Aus obiger Betrachtung folgt S=X i λiDi. (4.33) Gl. (4.33) kann wieder in Form der elastischen Steifigkeit ausgedrückt werden zu C=X i λ−1 iDi. (4.34) 4.2 Viskoelastische Materialmodelle für die Simulation akustischer Wellen in Polymeren 47 Spektrale Zerlegung für hexagonale Materialsymmetrie Es wird ein ideal elastisches und transversal isotropes Material angenommen. Die Extrusionsrichtung wird im Folgenden als Vorzugsrichtung (Index L) und die dazu orthogonale Fläche als IsotropieEbene (Index T) betrachtet. Im kartesischen (später auch im zylindrischen) Koordinatensystem wird die Vorzugsrichtung in z -Richtung angesetzt. Die elastische Nachgiebigkeitsmatrix S wird mittels fünf linear unabhängiger Konstanten beschrieben zu S=  1/ET−νT/ET−νL/EL000 1/ET−νL/EL000 1/EL000 1/µL0 0 sym. 1/µL0 1/µT  ;µT=ET 2(1+νT). (4.35) Entsprechend der elastischen Forderung gelte ET , EL , νT , νL , µL∈R+ . Die Zerlegung der Materialmatrix ergibt vier unterschiedliche Eigenwerte [The00; TP90] λ1=1 µT;λ2=1 µL λ3,4 =1−νT 2ET+1 2EL±v t"1−νT 2ET−1 2EL#2 +2ν2 L E2 L , (4.36) wobei λ1 und λ2 jeweils zweifache Eigenwerte darstellen. Die zugehörigen Projektionsmatrizen ergeben sich zu D1=  1−10000 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1  (4.37) D2=  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0  (4.38) 48 4 Wellenausbreitung in viskoelastischen Medien D3=  1 2cos2ωe1 2cos2ωe1 √2sin ωecos ωe0 0 0 1 2cos2ωe1 2cos2ωe1 √2sin ωecos ωe0 0 0 1 √2sin ωecos ωe1 √2sin ωecos ωesin2ωe0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  (4.39) D4=  1 2sin2ωe1 2sin2ωe−1 √2sin ωecos ωe0 0 0 1 2sin2ωe1 2sin2ωe−1 √2sin ωecos ωe0 0 0 −1 √2sin ωecos ωe−1 √2sin ωecos ωecos2ωe0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  , (4.40) wobei der Eigen-Winkel ωedefiniert wird zu tan 2ωe=−2√2νL (1−νT)EL ET−1. (4.41) Der Eigen-Winkel ist definiert im Intervall ωe∈[ 0, 180 ◦]⊂R und kann genutzt werden, um die Anisotropie des transversal isotropen Körpers zu charakterisieren. Für ωe∈[ 0 ◦ , 90 ◦] lässt sich zeigen, dass νL≤ 0 gilt [The00], was für die betrachteten Materialklassen i.A. nicht gilt und im Folgenden ausgeschlossen werden soll. Entsprechend wird ωe auf das Intervall [ 90 ◦ , 180 ◦] begrenzt, welches sich wiederum in drei Gebiete unterteilen lässt. Das Segment ωe∈[ 90 ◦ , 125,27 ◦] ist charakterisiert durch eine schwache Symmetrieachse σz bezogen auf die isotrope Ebene σx− σy (weak-axis material). Hingegen zeigt sich für das Segment ωe∈[ 135 ◦ , 180 ◦] eine starke Symmetrieachse (strong-axis material). Der Zwischenbereich mit ωe∈[ 125,27 ◦ , 135 ◦] ist durch die (im Intervall näherungsweise und an den Intervallgrenzen exakte) Entkopplung der Eigenwerte λ3 und λ4 hinsichtlich Longitudinal-Scherbewegung und Longitudinal-Kompressionsbewegung charakterisiert. Wird Gl. (4.41) an ωe= 125,27 ◦ ausgewertet, zeigt sich die Eigenschaft des sogenannten quasi-isotropen Verhaltens mit 1−νL EL=1−νT ET. (4.42) Die Eigenwerte und die Projektionsmatrizen vereinfachen sich dann zu λ3,qi =1+νL EL;λ4,qi =1−2νL EL=1 3KL(4.43) 4.2 Viskoelastische Materialmodelle für die Simulation akustischer Wellen in Polymeren 49 mit KLals longitudinaler Kompressionsmodul und D3,qi =1 6  1 1 −2 0 0 0 1 1 −2 0 0 0 −2−2 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  (4.44) D4,qi =1 3  1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  . (4.45) Ein Spezialfall für diesen Eigen-Winkel ist die Isotropie mit EL=ET , νL=νT und µL=µT . Für den anderen Grenzwert des Segments mit ωe= 135 ◦ muss der Nenner von Gl. (4.41) Null ergeben. Nach einigen elementaren Umformungen erhält man für die zugehörigen Eigenwerte 3 und 4 λ3;135◦=1 EL;λ4;135◦=1 EL. (4.46) Entsprechend ist auch für diesen Fall eine Entkopplung von Kompressionsund Scherkomponenten gegeben. Neben dem Eigen-Winkel ωe kann ebenso µL zur Charakterisierung der Anisotropie des Körpers dienen. Es existiert entsprechend eine Korrelation zwischen den beiden Größen, welche Theocaris phänomenologisch zu EL 2µL=aemp (1+|cot 2ωe|)=aemp 1+ 1 2EL−1−νT 2ET √2νLE−1 T(4.47) bestimmt hat [TP90; The93; The00], wobei er auf Basis der ihm zur Verfügung stehenden Messdaten, siehe Abb. 4.4(a), die Konstante zu aemp =√2 bestimmt hat. Durch Wahl dieser Konstanten tritt das Problem auf, dass bei Annahme eines isotropen Körpers mit EL=ET sowie νL=νT die Schätzung für µL nicht µT entspricht, also kein isotropes Material geschätzt wird. Durch Wahl der Konstanten zu aemp,kor =2√2 1+2√2(1+νL)(4.48) kann dies korrigiert werden [Rau12]. In Abb. 4.4 ist das Verhältnis EL 2µL über den Eigen-Winkel ωe aufgetragen. Die rote unterbrochene 50 4 Wellenausbreitung in viskoelastischen Medien Kurve beschreibt die Schätzung nach Theocaris mit aemp =√2 , die blaue unterbrochene Kurve die vorgeschlagene korrigierte Schätzung mit νL= 0,1 und die grüne unterbrochene Kurve die vorgeschlagene korrigierte Schätzung mit νL= 0,5. In Abb. 4.4(a) sind verschiedene Materialgruppen dargestellt, welche hexagonale Symmetrie aufweisen. Als rote Sterne sind (uniaxiale) Carbonund Graphit-Komposite eingezeichnet, welche sich durch eine starke Vorzugsrichtung auszeichnen. Als grüne Kreise sind gewebte Faser-Komposite eingezeichnet, welche eine deutlich weniger starke Anisotropie aufweisen. Ebenfalls dargestellt in schwarzen Rauten sind anorganische Kristalle der hexagonalen Klasse, welche sich nahe der quasi-isotropen Achse verorten lassen. Die genutzten Materialparameter sind [The98; The00] entnommen. In Abb. 4.4(b) sind hingegen die Eigenschaften verschiedener Polymere dargestellt. Die mit Sternen dargestellten Einträge sind nichtoder nur wenig gestreckte Polymere, die unterschiedlichen Farben lassen sich unterschiedlichen Literaturquellen zuordnen, siehe Anhang B.1. Die in Kreisen dargestellten Einträge bezeichnen stark gestreckte Polymere, sodass hier eine deutlich stärker ausgeprägte Vorzugsrichtung ensteht. 100 120 140 160 180 0 5 10 15 20 25 30 Eigen-Winkel ωe/◦ EL/2µL weak-axis strong-axis ←Quasi-Isotropie (a) Gewebte Faser-Komposite (o), Carbonund Graphit-Komposite (*) sowie anorganische (hexagonale) Kristalle (). 100 120 140 160 180 0 5 10 15 20 25 30 Eigen-Winkel ωe/◦ EL/2µL weak-axis strong-axis ←Quasi-Isotropie (b) Gestreckte (o) und (nahezu) nicht gestreckte Polymere (*). Abbildung 4.4: Darstellung des Verhältnisses EL 2µL über den Eigen-Winkel ωe . Unterbrochene Linien zeigen Schätzungen nach Gl. (4.47) mit aemp =√2 (rot) sowie mit aemp,kor für νL= 0,1 (blau) und νL=0,5 (grün). Es zeigt sich in beiden Grafiken, dass die Näherungsgleichung über viele Materialgruppen hinweg gute Näherungen liefert. Es ist zu beobachten, dass nahezu alle Daten im Bereich zwischen den berechneten Schätzungen mit aemp und aemp,kor liegen. Des Weiteren ist zu erkennen, dass solange keine gestreckten Polymere untersucht werden, keine ausgeprägten Anisotropien zu erwarten sind, sondern mit einer recht schwach ausgeprägten Form der transversalen Isotropie gearbeitet werden kann. Der Wertebereich der Eigen-Winkel der dargestellten Polymere liegt im mittleren Segment mit ωe∈[ 125,27 ◦ , 135 ◦] . Dies lässt im Folgenden zu, von entkoppelten Kompressions- 4.2 Viskoelastische Materialmodelle für die Simulation akustischer Wellen in Polymeren 51 und Scherkomponenten auszugehen. Außerdem wird auf Grundlage dieses Vorwissens für die Schätzung von µLdie korrigierte empirische Konstante genutzt. Die diskutieren Eigenschaften sind im Folgenden wichtig für die Einschätzung von berechneten Startwerten für den inversen Ansatz sowie für eine begründete Reduzierung der unabhängigen Parameter des viskoelastischen transversal isotropen Materialansatzes, siehe Kapitel 5.6. Spektrale Zerlegung für isotrope Materialien Der Spezialfall der Isotropie soll an dieser Stelle ebenfalls kurz betrachtet werden. Die elastische Steifigkeitsmatrix sei beschrieben durch den Kompressionsund den Schermodul, vgl. Gl. (4.1), mit C=  K+4 3µK−2 3µK−2 3µ0 0 0 K+4 3µK−2 3µ0 0 0 K+4 3µ0 0 0 µ0 0 sym. µ0 µ  . (4.49) Die spektrale Zerlegung nach Theocaris zeigt, wie für diese Materialsymmetrie-Gruppe zu erwarten, drei Eigenwerte mit zugehörigen Projektionsmatrizen: D1=1 3  1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  ;λ−1 1=3K(4.50) D2=1 3  2−1−1 0 0 0 −1 2 −1 0 0 0 −1−1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  ;λ−1 2=2µ(4.51) 52 4 Wellenausbreitung in viskoelastischen Medien D3=  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1  ;λ−1 3=µ(4.52) Es ist offensichtlich, dass hier ebenfalls eine Entkopplung der Kompressionsund Scherkomponenten erfolgt, wobei hinsichtlich der Scherkomponenten zwischen reiner Scherbewegung ( λ3 ) und der verkoppelten Scherung bei Volumenbewegungen ( λ2 ) unterschieden wird. Die Eigenwerte zwei und drei können mit ihren Projektionsmatrizen entsprechend zusammengefasst werden zu C=3KD1+µ(2D2+D3). (4.53) Generalisierung auf viskoelastisches Verhalten Für den Übergang zu viskoelastischem Verhalten (hier nun direkt im Fourier-Bereich) werden nun die reell-wertigen Eigenwerte in komplexe Funktionen in ω überführt, wobei die eindimensionalen Modellansätze aus Kapitel 2.3.1 angesetzt werden können. Dadurch wird jeder entkoppelten Eigenbewegung des Materials ein individueller Verlustfaktor zugeschrieben. In Form der (viskoelastischen) Steifigkeit geschrieben, ergibt sich ˜ C(iω) = X i λ−1 i(iω)Di. (4.54) Wird die Struktur des fraktionalen Zener-Modells nach Gl. (2.46) angesetzt (analog kann auch jedes andere beliebige Modell angesetzt werden), ergibt sich ˜ C(iω) = X i λ−1 i 1+(iωτεi)νi 1+(iωτσi)νiDi. (4.55) Dieser Ansatz ist im Sinne einer energetischen Betrachtung zu begründen. Es wird die potentielle Energie Epot eines Körpers betrachtet, welche sich für den elastischen Fall ausdrücken lässt zu [TP90] Epot =1 2σTε=1 2X i λiσT iσi. (4.56) Gl. (4.56) beschreibt eine Zerlegung der gesamten gespeicherten potentiellen Energie in orthogonale Teil-Energiespeicher, welche mit den orthogonalen Eigen-Spannungszuständen des Körpers korrespondieren. Jedem Teil-Energiespeicher wird durch den in Gl. (4.54) beschriebenen Ansatz 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter 53 Verluste zugeordnet, sodass jede mögliche Form der Energiespeicherung individuell verlustbehaftet ist. Im Kontext der Arbeit von Douven et al. [DSJ89] fällt auf, dass trotz der fünf linear unabhängigen elastischen Parameter des transversal isotropen Materials nur vier Parameter (Eigenwerte) als verlustbehaftet betrachtet werden müssen 6 . Der fünfte Parameter (Eigen-Winkel) hat keinen Einfluss auf die Energiespeicher, sondern reguliert die räumliche Gestalt der Eigenbewegungen. Im Kontext der Arbeit von Carcione [Car90] wurde hier eine generalisierte Variante der Matrix-Zerlegung gefunden. 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter Die bisherigen Ausführungen zur Wellenausbreitung in viskoelastischen Materialien haben sich auf Wellenausbreitung unter Freiraumbedingungen bezogen. Für den konzeptionellen Ansatz der Untersuchung von geometrisch limitierten Materialproben in Gestalt eines axialsymmetrischen Wellenleiters, siehe dazu Kapitel 5.1, sind die zuvor aufbereiteten Materialmodelle in den Kontext geführter Wellenausbreitung in hohlzylindrischen längshomogenen Wellenleitern zu setzen. Im Folgenden wird zur Berechnung der Dispersionsdiagramme die SBFEM 7 genutzt. Die Ergebnisse werden im Anschluss mit der Methode der Modalen Expansion weiterverarbeitet, um die Anregungsund Empfangscharakteristik der eingesetzten Schallwandler zu modellieren. Zunächst wird von kartesischen Koordinaten (x , y , z) zu Zylinderkoordinaten (r , θ , z) übergegangen, vgl. Abb. 4.5(a). Die Spannungsund Deformationsvektoren nach Gl. (2.8) bzw. Gl. (2.6) schreiben sich entsprechend σ=[σrr,σθθ,σzz,σθz,σrz,σθr]T ε=[εrr,εθθ,εzz, 2εθz, 2εrz, 2εθr]T. (4.57) Aufgrund der räumlichen Anregung der Probe, vgl. Abb. 5.2, wird die Betrachtung auf longitudinale Wellenleitermoden beschränkt. Diese Beschränkung des methodischen Ansatzes zur Charakterisierung der Materialparameter einer Probe impliziert formal zunächst keine Restriktionen bezüglich des transversal isotropen Materialansatzes. Durch die Axialsymmetrie der Moden entfallen alle parti6 Dieses Ergebnis darf nicht allgemein auf alle Materialsymmetrien übertragen werden. So ist für den Fall eines isotropen Materials der Übergang von Materialkonstanten zu zeitabhängigen Funktionen für die zwei unabhängigen Parameter identisch mit dem Ergebnis der hier dargestellten Matrix-Zerlegung. 7 Der Kern der hier verwendeten SBFEM ist in einer Kooperation mit der Bundesanstalt für Materialforschung und -prüfung entstanden. Insbesondere ist hier Dr. Hauke Gravenkamp zu nennen. 60 4 Wellenausbreitung in viskoelastischen Medien Spannungskomponente σ(j) zz (r , ω) und Tangential-Spannungskomponente σ(j) rz (r , ω) charakterisiert. Es sei an dieser Stelle bemerkt, dass Fama [Fam72] und Fraser [Fra75] die Orthogonalität und Vollständigkeit der Wellenfelder der Moden im zylindrischen Wellenleiter mit unterschiedlichen Randbedingungen für Spannungen und Verschiebungen untersucht haben. Für die hier gewählten Randbedingungen konnte keine Vollständigkeit nachgewiesen werden, jedoch konnten Puckett und Petersen [PP05] zeigen, dass zumindest im Sinne kleinster Fehlerquadrate eine eindeutige Lösung erzielt werden kann. Im Folgenden wird sich daher an der Vorgehensweise von Puckett und Petersen orientiert. Empfängerseitig ist zu beachten, dass nur Normalspannungen zum elektrischen Empfangssignal beitragen. Der Empfangs-Schallwandler empfängt dabei integral über die Oberfläche der Stirnseite der Probe. Diese integrale Eigenschaft kann effizient in der Formulierung der Randbedingungen durch Normierung der modalen Spannungskomponenten σ(j) rz (r , ω) und σ(j) zz (r , ω) auf ihr jeweiliges Integral über den Wellenleiterquerschnitt berücksichtigt werden [PP05]. Entsprechend werden normierte Größen σ(j) zz,rz(r,ω) = σ(j) zz,rz(r,ω) 2πrri+Dw riσ(j) zz,rz(r,ω)rdr (4.79) für die Modale Expansion verwendet. In der räumlich und zeitlich diskreten Analyse müssen die aufgestellten Randbedingungen an diskreten Kreisfrequenzen ωn ausgewertet werden. An einer diskreten Kreisfrequenz existieren dann Nj Wellenleitermoden, deren normierte modale Normalund Tangential-Spannungskomponenten durch die radiale Diskretisierung räumlich diskret vorliegen. Die Diskretisierung ist abhängig von der Frequenz und in der Regel nicht äquidistant gewählt, sodass die Spannungskomponenten entlang des Radius auf Nr äquidistante Punkte interpoliert werden. Die Randbedingungen aus Gl. (4.77) und Gl. (4.78) können dann diskret formuliert werden zu  σz[1] σz[2] . . . σz[Nr] 0 0 . . . 0  | {z } V =  σ(1) zz [1,ωn]σ(2) zz [1,ωn]··· σ(Nj) zz [1,ωn] σ(1) zz [2,ωn]σ(2) zz [2,ωn]··· σ(Nj) zz [2,ωn] . . .. . .. . . σ(1) zz [Nr,ωn]σ(2) zz [Nr,ωn]··· σ(Nj) zz [Nr,ωn] σ(1) rz [1,ωn]σ(2) rz [1,ωn]··· σ(Nj) rz [1,ωn] σ(1) rz [2,ωn]σ(2) rz [2,ωn]··· σ(Nj) rz [2,ωn] . . .. . .. . . σ(1) rz [Nr,ωn]σ(2) rz [Nr,ωn]··· σ(Nj) rz [Nr,ωn]  | {z } Σωn · A(1)[ωn] A(2)[ωn] . . . A(Nj)[ωn]  | {z } A . (4.80) 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter 61 Zur Lösung wird das Gleichungssystem linksseitig mit ΣT ωn multipliziert und dann die Inverse von ΣT ωnΣωndurch QR-Zerlegung berechnet. Es ergibt sich A=ΣT ωnΣωn−1ΣT ωnV. (4.81) Die Performanz dieses Ansatzes der Modalen Expansion ist abhängig von der Anzahl der Moden Nj , welche zur Rekonstruktion der Randbedingung genutzt werden. Eine höhere Anzahl von Moden erhöht die Stabilität der Methode. Die geometrischen Abmessungen sowie die Beschränkung auf Polymere mit niedrigen Wellenausbreitungsgeschwindigkeiten (in der Folge kleine Wellenlängen), verglichen mit z.B. Metallen, gewährleistet eine ausreichende Anzahl von Moden für diesen Ansatz. Zur Abschätzung der Anzahl ausbreitungsfähiger Moden genügt die Betrachtung das Verhältnis von Mantelstärke des Wellenleiters zur halben Wellenlänge der Longitudinalwelle. Ist nun bekannt, mit welcher Amplitude eine Mode j an der Frequenz ωn angeregt wird, können die Moden durch den Wellenleiter von Position z= 0 an die Position z=l mithilfe der zugehörigen bekannten Wellenzahl kz( i ωn) propagiert werden. Es ergibt sich der Frequenzgang des Wellenleiters, hier als kontinuierliche Funktion in ωnotiert, zu HWL(iω) = X j A(j)(ω)exp (−ikz(iω)l). (4.82) 4.3.3 Validierung der Methode Zur Validierung des aufgestellten Modells werden die berechneten transienten Signale mit Ergebnissen aus transienten FEM-Simulationen verglichen. Es wird das für die Simulation von mechatronischen Sensoren und Aktoren ausgelegte Softwarepaket CAPA (WisSoft) [Wis05] genutzt. Da CAPA nur die klassische Rayleigh-Dämpfung als Dämpfungsmodell unterstützt, muss sich in diesem Abschnitt für die Validierung auf eben dieses Modell beschränkt werden. Eine Skizze des genutzten FEM-Modells ist in Abb. 4.8 dargestellt. Es werden lineare mechanische 2D axialsymmetrische Elemente genutzt, um die axialsymmetrische Struktur (Rotationsachse z ) räumlich zu diskretisieren. Die axiale Ausdehnung des Bilanzraums ist insgesamt drei Mal die Länge l des Wellenleiters, wobei der mittlere Part den Wellenleiter beschreibt (weiß) und die beiden Randbereiche (grau) Dämpfungsbereiche darstellen. In den Dämpfungsbereichen werden die Rayleigh-Dämpfungsparameter αM und αK logarithmisch erhöht, siehe ebenfalls Abb. 4.8. Diese Vorgehensweise bezüglich der Dämpfungsregionen bewirkt akkurate Ergebnisse bzw. vernachlässigbar kleine Reflexionen aus den Dämpfungsbereichen [Lan06]. An z= 0 werden loads bezüglich der Normal-Spannungskomponenten definiert und an z=l das Flächenintegral bezüglich der Normal-Spannungen ausgewertet. Die räumliche Diskretisierung in radialer Richtung beträgt 33 Knoten pro Wellenlänge bezüglich der longitudinalen Freiraum-Wellenausbreitungsgeschwindigkeit 62 4 Wellenausbreitung in viskoelastischen Medien bei 1,6 MHz . Die zeitliche Diskretisierung für den Newmark Zeitschritt-Algorithmus beträgt 12,5 ns . 100 10 1 αM,K(z) αM,K(z=0) r ri ri+DW z z=0 z=l Abbildung 4.8: Skizze des FEM Simulationsmodells zur transienten Simulation der Wellenausbreitung im axialsymmetrischen Wellenleiter. Betrachtungen auf Basis des komplexen Materialmodells: Für einen Vergleich zwischen der halb-analytischen Simulation (SBFEM und Modale Expansion) und einer numerischen Zeitschritt-Simulation in der FEM ist bezüglich der halb-analytischen Simulation ein komplexes Materialmodell nach Gl. (4.55) basierend auf dem Anti-Zener-Modell zu verwenden. Die Äquivalenz zwischen Rayleigh-Dämpfung und Anti-Zener Modell wurde bereits für den eindimensionalen Fall in Kapitel 4.2.1 diskutiert. Für den dreidimensionalen Fall ist jedoch eine Anpassung der Absorptionseigenschaften der Eigenbewegungen vorzunehmen. Nach Gl. (4.11) manifestiert sich die Absorption durch das komplexe Materialmodell (in der Einheit Np m−1 ) in Abhängigkeit der Wellenausbreitungsgeschwindigkeit. Bezogen auf die Eigenbewegung mit dem Eigenwert λj gilt für das Anti-Zener-Modell: αj(ω) = −ω1 cj Im       sω2 α2 M,j+ω2(1+i tan δj)−1       (4.83) tan δj=αK,jω+αM,j(αM,jαM,j+1) ω;cj=qλjρ−1(4.84) Das bedeutet, dass, obwohl verschiedenen Eigenbewegungen identische Verlustwinkel zugewiesen werden, die Absorption unterschiedlich ist. Das Rayleigh-Dämpfungsmodell wirkt im Zeitschrittverfahren derart, dass ein Abklingen über der Zeit für alle Wellenarten identisch stattfindet. Entsprechend sind die Parametrisierungen der Anti-Zener-Modelle für jede Eigenbewegung so anzupassen, dass effektiv für alle Eigenbewegungen 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter 63 eine identische Laufzeit-basierte Absorption ˜αj(ω) = αj(ω)cj vorliegt. Dazu wird ein einfacher inverser Ansatz gewählt, in welchem die Parameter αM,j,k und αK,j,k durch Optimierung mit dem Nelder-Mead Simplex-Algorithmus so variiert werden, dass der Verlauf der Absorption ˜αj=k(ω) hinreichend gut approximiert wird. Es existiert dazu keine geschlossene Lösung, da die Skalierung der Absorption frequenzunabhängig durchzuführen wäre, dies jedoch mit dem Modellansatz nicht realisiert werden kann. Entsprechend sind bei diesem Ansatz der komplexen Modellierung für den Vergleich mit den FEM-Simulationen im Zeitschrittverfahren Abweichungen zu erwarten. Das vorgeschlagene Vorgehen resultiert jedoch in hinreichenden Ergebnissen für eine Überprüfung des komplexen Materialmodell-Ansatzes. Ein alternativer Ansatz, welcher sich näher an der Natur der Rayleigh-Dämpfung orientiert, wird im nächsten Abschnitt betrachtet. Als Testszenario wird ein Wellenleiter mit Innenradius ri= 3 mm und einer Mantelstärke von Dw= 6,5 mm gewählt. Simuliert wird die Propagation einer Gaußpuls-förmigen Anregung mit einer Mittenfrequenz von 1 MHz und ca. 60% relativer Bandbreite. Es werden vier unterschiedliche Szenarien betrachtet: •Ideal elastische Eigenschaften bei isotroper Materialsymmetrie, siehe Abb. 4.9 •Viskoelastische Eigenschaften bei isotroper Materialsymmetrie, siehe Abb. 4.10 •Ideal elastische Eigenschaften bei hexagonaler Materialsymmetrie, siehe Abb. 4.11(a) •Viskoelastische Eigenschaften bei hexagonaler Materialsymmetrie, siehe Abb. 4.11(b) Die gewählten Materialparameter für den isotropen Fall sind in Tabelle 4.2 aufgeführt. Die MaterialKennzeichnungen ID1 und ID2 sind an Materialparametern von PEEK und PP angelehnt. Der transversal isotrope elastische Materialparametersatz in Tabelle 4.3 beschreibt PE Methacrylat in Anlehnung an den Datensatz aus [Hel79]. Tabelle 4.2: Tabelle isotroper Materialparameter für den Vergleich mit FEM-Simulationen. Material K/GPa µ/GPa ρ/kg m−3αM/sαK/s−1 ID1 (PEEK) 6,8369 1,7193 1300 11 ·1042,8 ·10−9 ID2 (PP) 4,4575 1,5738 906 40 ·1045·10−9 Tabelle 4.3: Tabelle transversal isotroper Materialparameter mit quasi-isotroper Annahme für den Vergleich mit FEM-Simulationen. Material EL/GPa ET/GPa µL/GPa νLρ/kg m−3αM/sαK/s−1 ID3 (PE) 7,14 5,41 2,29 0,31 950 15 ·1048·10−9 64 4 Wellenausbreitung in viskoelastischen Medien Die Ergebnisse für die ideal elastischen isotropen Materialparameter sind in Abb. 4.9 zusammengestellt. Neben dem Ergebnis der FEM-Simulation (rot) und der SBFEM-Simulation (blau) ist zusätzlich der Eintreffzeitpunkt der Rayleigh-Welle (grün) eingezeichnet. Die Oberflächenwelle wird in der halb-analytischen Simulation nicht berücksichtigt, was zu Abweichungen ab diesem Zeitpunkt zwischen den Simulationen führt. Des Weiteren sind bei der halb-analytischen Simulation sehr kleine Signalanteile im Bereich bis 5 µ s zu sehen, welche durch die zyklische Rücktransformation von der Frequenzbereichs-Modellierung in die Zeitdarstellung zu erklären sind. Sie beschreiben Signalanteile, welche formal nach Ende der Simulationszeit Tsim =∆f−1 sim = 100 µ s eintreffen würden (kleine Gruppengeschwindigkeit) und durch die zyklische FFT an den Anfang das Empfangssignals projiziert werden. Diese Signalanteile werden nur in diesem Abschnitt zur Diskussion dargestellt, in der folgende Arbeit jedoch ausgeblendet. Abgesehen von den genannten Abweichungen ist die Übereinstimmung der simulierten Signale gut. Dies gilt insbesondere vor dem Hintergrund, dass die Rechenzeit 9 der halb-analytischen Simulation für das reelle Material ID1 20,5 s und für das reelle Material ID2 17,5 s 10 beträgt. Im Vergleich dazu wird für die FEM-Simulation ca. 30 Minuten benötigt 11 . Der benötigte Arbeitsspeicher beträgt ca. 400 MB . Nur durch die beschleunigte Simulation wird das in dieser Arbeit vorgeschlagene inverse Verfahren zur Materialparameterbestimmung anwendbar. Durch die Erhöhung der betrachteten Signalfrequenzen steigt die Anzahl der Freiheitsgrade in der FEM weiter an und führt zu noch größeren Simulationszeiten. Bspw. dauert eine FEM-Simulation einer Transmission eines Pulses mit 2 MHz Mittenfrequenz durch den oben beschriebenen Wellenleiter bereits deutlich länger als eine Stunde bei einem Speicherbedarf von ca. 2700 MB. Zeit t/µs Normiertes Signal 0 5 10 15 20 -1 -0.5 0 0.5 1 (a) ID1; Isotrop; Elastisch. Zeit t/µs Normiertes Signal 0 5 10 15 -1 -0.5 0 0.5 1 (b) ID2; Isotrop; Elastisch. Abbildung 4.9: Vergleich zwischen transienter FEM-Simulation (rot) und SBFEM-Simulation (blau). Eintreffzeitpunkt der Rayleigh-Welle (grün). 9 Simulationen werden ausgeführt auf einem Standard PC mit Intel(R) Core(TM) i7-4770K CPU mit 3,5 MHz Taktrate und 16 GB RAM. 10Frequenzbereich der Simulation: 0 Hz −2,5 MHz in 1/100 MHz Schritten. Zunächst keine Parallelisierung. 11Mittlerer Wert für die Rechenzeit der verschiedenen FEM-Simulationen 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter 65 Wird nun die Rayleigh-Dämpfung in der FEM-Simulation bzw. das Anti-Zener-Modell mit der Rayleigh-Dämpfungs-Interpretation für die SBFEM-Simulation betrachtet, ergeben sich die in Abb. 4.10 zusammengestellten Ergebnisse. Die Simulationsergebnisse passen weiterhin gut zusammen. Es ist allerdings auffällig, dass die nicht-kausalen Signalanteile etwas zunehmen. Dies ist durch die im Anti-Zener-Modell berücksichtigte Materialdispersion zu erklären. Hin zu kleinen Frequenzen sinkt die simulierte Phasengeschwindigkeit der Freiraum-Welle sehr stark ab, was zu langsamen Wellenanteilen im Empfangssignal führt, welche wiederum durch die zyklische FFT an den Anfang des simulierten Signals projiziert werden. Die durch die Materialdispersion verursachte Verringerung der Phasengeschwindigkeit findet -aufgrund der Parametrisierung der Modelle für das Material ID2schon bei höheren Frequenzen statt als bei dem Material ID1. Dies spiegelt sich in den Zeitsignalen durch die Ausprägung der entsprechenden langsamen Signalanteile wieder. Die Rechenzeiten unter Zugrundelegung der komplex-wertigen Materialien betragen 29 s für das Material ID1 und 21 s für das komplexe Material ID2. Zeit t/µs Normiertes Signal 0 5 10 15 20 -1 -0.5 0 0.5 1 (a) ID1; Isotrop; Rayleigh-Dämpfung. Zeit t/µs Normiertes Signal 0 5 10 15 -1 -0.5 0 0.5 1 (b) ID2; Isotrop; Rayleigh-Dämpfung. Abbildung 4.10: Vergleich zwischen transienter FEM-Simulation (rot) und SBFEM-Simulation (blau). Abb. 4.11 zeigt die elastische und die viskoelastische halb-analytische Simulation des transversal isotropen Materials ID3 jeweils im Vergleich zur FEM-Simulation. Da die Rechenzeit im Wesentlichen von der Anzahl der berechneten Moden abhängt, erhöht sich diese nicht durch die verallgemeinerte Materialsymmetrie und beträgt für den elastischen Fall ca. 16 s und für den viskoelastischen Fall ca. 19 s. Die Übereinstimmung der simulierten Signale ist auch hier hinreichend gegeben. Erweitere Modellierung hinsichtlich Rayleigh-Dämpfungsmodell-Äquivalenz: Wie im vorherigen Absatz beschrieben, ist ein Vergleich des komplex-wertigen Materialmodells in der SBFEM mit dem Zeitschrittverfahren in der FEM nur eingeschränkt möglich. Durch die Geometrie des Wellenleiters sowie dessen räumliche Anregung an der Stirnfläche mit Normalspannungen lässt sich ein entkoppelter Ansatz von Materialdispersion im Sinne der frequenzabhängigen Wellenaus- 66 4 Wellenausbreitung in viskoelastischen Medien Zeit t/µs Normiertes Signal 0 5 10 15 -1 -0.5 0 0.5 1 (a) ID3; Transversal isotrop; Elastisch. Zeit t/µs Normiertes Signal 0 5 10 15 20 -1 -0.5 0 0.5 1 (b) ID3; Transversal isotrop; Rayleigh-Dämpfung. Abbildung 4.11: Vergleich zwischen transienter FEM-Simulation (rot) und SBFEM-Simulation (blau). breitungsgeschwindigkeit, und der frequenzabhängigen Absorption formulieren, vgl. [BGR15b]. In dieser entkoppelten Variante kann die frequenzabhängige Absorption im Sinne der RayleighDämpfung direkt in Form der Absorption pro Zeiteinheit angesetzt werden. Zunächst wird eine reell-wertige Formulierung der Materialmatrix gesucht, welche dieselbe Frequenzabhängigkeit aufweist wie eine komplex-wertige Materialmatrix auf Basis des RayleighDämpfungsmodells nach Gl. (4.21). Ausgehend von der komplexen Formulierung ˜ C(iω) = 1 1−α2 M (iω)2 1+iω −αM (iω)2(1+αMαK)+αK!!C(4.85) =γAZ(iω,αM,αK)C, (4.86) ergibt sich bspw. die Longitudinalwellengeschwindigkeit in z-Richtung cph,L,z(ω) = qC33 ρ−1Re        qγAZ(iω,αM,αK)!−1       −1 , (4.87) wobei C33 ∈R+ das 33-Element aus Matrix C beschreibt. Der Ausdruck für cph,L,z(ω) ist eine reelle Funktion der Frequenz, sodass eine äquivalente reell-wertige Formulierung angesetzt werden kann, um dieselbe Frequenzabhängigkeit zu beschreiben. Dies ist für ˘ C33(ω) = C33 Re        qγAZ(iω,αM,αK)!−1       −2 (4.88) 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter 67 erfüllt. Diese Herleitung gilt für jedes Element aus Matrix C, sodass geschrieben werden kann ˘ C(ω) = CRe        qγAZ(iω,αM,αK)!−1       −2 . (4.89) Zugehörig zur reell-wertigen frequenzabhängigen Materialmatrix ˘ C(ω) steht die Absorption des Materials. Diese ist im Freiraum nach Gl. (4.11) beschrieben durch αC33 (ω) = −ω1 pC33 ρ−1Im        qγAZ(iω,αM,αK)!−1       . (4.90) Durch Projektion der streckenabhängigen Absorption αC33 (ω) auf die zeitabhängige Absorption ˜α(ω)entfällt die Abhängigkeit der Absorption von der Steifigkeit ˜α(ω) = αC33 (ω)cph,L,z(ω) =−ωRe        qγAZ(iω,αM,αK)!−1       −1 Im        qγAZ(iω,αM,αK)!−1       . (4.91) Es erfolgt zur Simulation in der SBFEM zunächst die Berechnung der Dispersionsdiagramme auf Basis der Materialmatrix aus Gl. (4.89) sowie die unveränderte Berechnung der modalen Amplituden nach Gl. (4.80). Bei der Formulierung des Frequenzgangs des Wellenleiters nach Gl. (4.82) muss nun nachträglich der Einfluss der Absorption wieder eingeführt werden. Die Impulsantwort des noch ungedämpften Wellenleiters als inverse Fourier-Transformierte des Frequenzgangs sei beschrieben als ˘ h(t) . Jede harmonische Komponente von ˘ h(t) muss nun nachträglich um den Faktor exp(−˜α(ω)t) gedämpft werden. Jede harmonische Schwingung von ˘ h(t) zeigt sich in der Fourier-Reihe ˘ h(t) = X ι Hιexp(−iωιt)mit ωι=ι2πT−1,ι∈N, (4.92) wobei Hι∈C die komplexen Fourier-Koeffizienten und T∈R+ die Periodizität von ˘ h(t) (zyklische diskrete Transformation) beschreibt. Werden nun die ungedämpften harmonischen Schwingungen durch gedämpfte Schwingungen ersetzt und die Gleichung zur Einhaltung der Kausalität mit der Sprungfunktion (Heavisidefunktion) H(t)multipliziert, ergibt sich [Rau12; BGR15b] hdamp(t) = X ι HιH(t)exp(−˜αιt)exp(−iωιt)mit ˜αι=˜α(ωι). (4.93) 68 4 Wellenausbreitung in viskoelastischen Medien Da die Modellierung des Frequenzgangs im Fourier-Bereich stattfindet, wird die Impulsantwort nach Gl. (4.93) wieder transformiert zu HWL,damp(iω) = X ιFH(t)exp(−˜αιt)∗[Hιδ(ωι) + H∗ ιδ(−ωι)]. (4.94) Wird die Faltungsoperation ausgeführt, ergibt sich HWL,damp(iω) = X ι Hι 1 ˜αι+i(ω−ωι)+H∗ ι 1 ˜αι+i(ω+ωι)!, (4.95) was nun den Frequenzgang unter Berücksichtigung der Absorptionseigenschaft des Wellenleitermaterials beschreibt. Durch diesen entkoppelten Ansatz von Absorption und frequenzabhängiger Steifigkeit werden implizit einige Approximationen eingeführt, welche im Folgenden diskutiert werden. • Die Berechnung der Dispersionseigenschaft des Wellenleiters wird allein auf Grundlage der reell-wertigen (aber frequenzabhängigen) Materialmatrix durchgeführt, vgl. Gl. (4.89). • Es wird nachträglich die Absorption des Materials unter Freiraumbedingungen in die Modellierung einbezogen, jedoch beinhaltet dies implizit eine Vernachlässigung der geometrischen Dämpfung der Wellenleiterstruktur. Ein Vergleich zwischen komplexer Modellierung und der entkoppelten Modellierung zeigt, dass die getroffenen Vereinfachungen eine gute Näherung darstellen, siehe Abb. 4.12(b-c). Die Phasengeschwindigkeiten der einzelnen Moden verändern sich nur wenig durch Einführung komplexer Materialparameter. Dies stimmt mit der Beobachtung für analytische komplex-wertige Wellenleitermodellierung überein, bei welchem das reell-wertige Problem zunächst genutzt wird, um gute Startwerte für die Suche nach der komplex-wertigen Wellenzahl zu erhalten (Nullstellensuche in der komplexen Ebene) [BBM14; XJQ05; PLA97]. Des Weiteren zeigt sich, dass für die vorliegende Wellenleitergeometrie in Kombination mit der Anregung von longitudinalen Wellenleitermoden über die Stirnfläche des Wellenleiters jene Wellenleitermoden angeregt werden, deren Phasengeschwindigkeit nahe der longitudinalen Freiraum-Wellengeschwindigkeit liegen, siehe Abb. 4.12(a). Ebenfalls ist zu erkennen, dass eben diese angeregten Bereiche der Wellenleitermoden eine minimale geometrische Dämpfung aufweisen und somit nahe an der Freiraumwellenabsorption liegen. Die Ergebnisse der entkoppelten Simulation und FEM-Simulationen sind in Abb. 4.13 für die Materialparameter ID1, ID2 und ID3 dargestellt. Die Rechenzeiten für die Simulationen sind vergleichbar mit den Simulationen der rein elastischen Materialmatrix. Die Nachbearbeitung des Frequenzgangs ist entsprechend deutlich weniger aufwendig als die Berechnung der Eigenwerte. Beim Vergleich der Ergebnisse wird deutlich, dass die entkoppelte Variante leichte Vorteile gegenüber der komplex- 4.3 Simulation transienter Signale im viskoelastischen Wellenleiter 69 Phasengeschwindigkeit c(j) ph /m s−1 Frequenz f/kHz 0 500 1000 1500 2000 2500 0 1000 2000 3000 4000 5000 6000 0 0.2 0.4 0.6 0.8 1 A(j) (a) Phasengeschwindigkeit nach komplexer Modellierung der Materialmatrix Phasengeschwindigkeit c(j) ph /m s−1 Frequenz f/kHz 0 500 1000 1500 2000 2500 0 1000 2000 3000 4000 5000 6000 0 0.2 0.4 0.6 0.8 1 A(j) (b) Vergleich zwischen entkoppelter (grau) und komplexer (farbig) Modellierung. Bezüglich der komplexen Rechnung sind nur die Moden dargestellt, deren Absorption kleiner als 0.541 Np/mm ist. Absorption α/Np/mm Frequenz f/kHz 0 500 1000 1500 2000 2500 0 0.05 0.1 0.15 0.2 0.25 0 0.2 0.4 0.6 0.8 1 A(j) (c) Absorption der Moden j. Abbildung 4.12: Vergleich zwischen entkoppelter Simulation und komplexer Modellierung der Materialmatrix. Farbkodierung der Kreise entspricht den modalen Amplituden A(j)(ω) . Die rot und grün unterbrochenen Linien zeigen die Freiraum-Longitudinalwellengeschwindigkeit und Transversalwellengeschwindigkeit sowie deren Freiraumwellenabsorption. 76 5 Realisierung des inversen Messverfahrens 5.2.1 Thermokammer und Temperaturregelung Die Thermokammer hat quadratische Außenabmessungen mit einer Seitenlänge von ca. 400 mm . Zur Isolation ist der Innenraum mit ca. 30 mm dicken Platten aus Polystyrol ausgelegt. An den Kälte-Umwälzthermostat ist ein Radiator (EK Water Blocks EK-CoolStream RAD XTX 240) angeschlossen, welcher mit zwei 120 mm Lüftern bestückt ist. Seitlich des Radiators gewährleisten aufrecht stehende Bleche eine Luftzirkulation innerhalb der Thermokammer. Der Pt100 Temperatursensor, welcher zur Temperaturregelung genutzt wird, ist an einem der Leitbleche so angebracht, dass der Sensor die Temperatur im Bereich der Messstrecke erfasst, vergleiche Abb. 5.5 (a). Der Kälte-Umwälzthermostat arbeitet mit deionisiertem Wasser als Temperierflüssigkeit, sodass der realisierbare Temperaturbereich laut Datenblatt 5 ◦C bis 90 ◦C beträgt. In Abb. 5.5(b) ist das Verhalten der Thermokammer aufgezeigt, welches sich als PT1-Verhalten gut beschreiben lässt. Die vorzunehmenden Messungen in der Thermokammer werden im Folgenden jeweils 30 Minuten nach Erreichen der Soll-Temperatur durchgeführt, um eine homogene Temperaturverteilung in der Probe und den Schallwandlern zu erreichen, siehe Anhang B.2. (a) Skizze der Thermokammer. Zeit t/min Temperatur T/◦C 0 100 200 300 20 25 30 35 40 45 50 (b) Temperaturtreppe; schwarz: Soll-Temperatur, blau: Ist-Temperatur. Abbildung 5.5: Thermokammer. 5.2.2 Handyscope HS5 Das Handyscope HS5 ist ein Zweikanal USB-Oszilloskop sowie ein Einkanal arbitrary waveform generator (AWG), siehe Abb. 5.6. Oszilloskop und AWG können über die USB-Schnittstelle aus MATLAB heraus angesteuert werden. Für die im Folgenden aufgeführten Messungen werden konstant folgende zentrale Einstellungen verwendet: 5.2 Konzeption und Funktionalität des Versuchsaufbaus 77 Abbildung 5.6: Handyscope HS5 vom Unternehmen Tiepie engineering, Niederlande [Tie14]. •Abtastrate: fs,ADU =fs,DAU =50 MHz •Quantisierung: NADU =NDAU =14 bit •Triggerereignis der Analog-Digital-Umsetzung: Generatorstart Zentrale Punkte zur Bewertung des Messsystems hinsichtlich präziser, reproduzierbarer Messungen sind der aperture jitter tj,aperture der clock jitter tj,clock sowie der trigger jitter tj,trigger . Seitens der Analog-Digital-Umsetzung sind zunächst die aperture jitter und der clock jitter von Bedeutung. Der clock jitter wird im Datenblatt über die Genauigkeit des Quarzoszillators mit ± 0.0001 % definiert. Bei einer Abtastrate von 50 MHz ergibt sich daraus ein clock jitter von tj,clock = 0,02 ps . Der i.A. deutlich dominantere aperture jitter wird im Datenblatt nicht angegeben. Beide Jitter, zum RMS-Wert tj,RMS zusammengefasst, stören die Annahme einer zeitlich äquidistanten diskreten Aufzeichnung. Dieses Phänomen kann als Phasenmodulation modelliert werden, wobei die Phase als stochastischer Prozess der Zeit angenommen wird [PMC11]. Die Phasenmodulation wird im quantisierten Signal erst dann sichtbar, wenn sie so groß ist, dass der quantisierte Signalwertunterschied zwischen phasenmoduliertem und originalem Signalwert gerade eine Quantisierungsstufe übersteigt. Damit lässt sich eine obere Grenze für den erlaubten Jitter in Abhängigkeit von der maximal erlaubten Signalsteigung (maximale Signalfrequenz, fmax) und der Anzahl der verwendeten Bits zur Quantisierung angeben [PMC11]: tj,RMSmax =1 2πfmax2NADU−1. (5.5) Mit NADU = 14 bit sowie einer Begrenzung der zu erwartenden maximalen Signalfrequenz auf fmax = 4 MHz ergibt sich ein maximal zulässiger Jitter von tj,RMSmax ≈ 4,86 ps . Aus dem RMSWert lässt sich wiederum eine zulässige Grenze für den unbekannten aperture jitter bestimmen, welcher, aufgrund des sehr kleinen clock jitter, mit dem maximalen RMS-Jitter näherungsweise übereinstimmt. Wird diese Grenze nicht eingehalten, zeigt sich in den digitalisierten Messdaten ein frequenzabhängiger Signal-zu-Störabstand (SNR) mit SNR =20log10 1 2πf tj,RMS !, (5.6) 78 5 Realisierung des inversen Messverfahrens welcher hinsichtlich einer Unsicherheitsanalyse im inversen Messverfahren nur schwer zu berücksichtigen ist. Aktuelle ADU weisen einen aperture jitter von zumeist < 1 ps auf [Kes08; LF03], sodass im Folgenden davon ausgegangen wird, dass es hinsichtlich der ADU-seitigen Umsetzung bei 50 MHz Abtastrate und 14 bit Quantisierung keine Jitter-Probleme zu beachten gilt. Der trigger jitter des Handyscope HS5 wird für den Fall der Triggerung auf den Generatorstart mit ≤ 4 Abtastwerten angegeben. Hinsichtlich einer Absolutlaufzeitmessung, sowie einer Mittelung über mehrere Empfangssignale im Zeitbereich, ist hier eine externe Stabilisierung des Systemverhaltens notwendig. Zunächst wird der trigger jitter experimentell für die oben aufgelisteten Einstellungen am Handyscope näher untersucht. Dazu werden Kanal 1 und Kanal 2 des Handyscope mit gleich langen BNC Leitungen mit dem AWG des Handyscopes verbunden. Beiden Kanälen wird ein Wertebereich von ± 8V zur Quantisierung zugewiesen. Dem Generator wird vorgegeben, einen gaußpulsmodulierten Sinus (im Folgenden kurz Gaußpuls) mit einer Mittenfrequenz fc= 1 MHz und relativer Bandbreite B= 110 % sowie einem Spitzenwert von 3,5 V auszugeben. Zunächst sei der Effekt des trigger jitter im Zeitbereich dargestellt, siehe Abb. 5.7. Es sind 50 Wiederholmessungen zu sehen, wobei leicht zu erkennen ist, dass die auf Kanal 1 aufgezeichneten Signalverläufe hinsichtlich der Zeitbasis schwanken. Zur näheren Analyse des Jitters wird die jeweilige Position des Signalmaximums für jede Einzelmessung bestimmt und die Stichprobe auf 10 4 Messungen erhöht. Zunächst wird die Position des Signalmaximums abtastwertgenau bestimmt und nachfolgend der Verlauf des Gaußpulses um das Maximum mit einem Polynom 2ten Grades approximiert, wobei jeweils ein Wert vor und ein Wert nach dem Maximalwert berücksichtigt wird. Basierend auf dem Polynomausdruck wird dann analytisch die Position des Maximums bestimmt, sodass mit Zwischenwertgenauigkeit eine Angabe über die Lage des Signalmaximums gemacht werden kann. Abb. 5.8(a) und (b) zeigen die relative Häufigkeitsdichte der Position des Signalmaximums jeweils bestimmt auf Kanal 1 und Kanal 2. Ebenfalls eingezeichnet ist die Position des Soll-Signalmaximums, wie es vom Generator vorgegeben ist (rote Linie). Es zeigt sich eine zufällige Verteilung der bestimmten Positionen in einem Bereich von ca. 4 − 5 Abtastwerten, wobei zu erkennen ist, dass der jeweilige Mittelwert der beiden Verteilungen etwas geringer ist als die Vorgabe. Das dargestellte Verhalten ist für eine präzise Laufzeitmessung nicht angemessen. Aus diesem Grund wird ein Kompensationsverfahren des trigger jitters angewendet, welches darauf beruht, dass der Jitter auf Kanal 1 und Kanal 2 korreliert ist. In Abb. 5.9(a) ist ein Ausschnitt aus den gemessenen Daten dargestellt. Auf der Ordinate ist die berechnete Position des Signalmaximums und auf der Abszisse der Zählindex der jeweiligen Messung dargestellt. Durch die Berechnung des Korrelationskoeffizienten zwischen den berechneten Positionen der Signalmaxima zwischen Kanal 1 und Kanal 2 zu 0,989 wird die Annahme der Korreliertheit bestätigt. Zur Kompensation des trigger jitters sowie der systematischen Abweichung von der Generatorvorgabe, wird das Generatorsignal auf Kanal 1 aufgezeichnet. Kanal 2 kann dann für die Datenerfassung der eigentlichen 5.2 Konzeption und Funktionalität des Versuchsaufbaus 79 Abtastwert n Spannung U/V 60 80 100 120 140 -1 0 1 2 3 Abbildung 5.7: 50 Messungen bei identischen Einstellungen des Handyscopes. Abtastwert des Signalmaximums nmax Rel. Häufigkeitsdichte 100 102 104 106 0 0.1 0.2 0.3 0.4 0.5 (a) Kanal 1. Abtastwert des Signalmaximums nmax Rel. Häufigkeitsdichte 100 102 104 106 0 0.1 0.2 0.3 0.4 0.5 (b) Kanal 2. Abbildung 5.8: Relative Häufigkeitsdichte der Position des Signalmaximums in Abtastwerten. Die rote Linie zeigt die am Generator vorgegebene Position des Signalmaximums an. Messinformation genutzt werden. Die auf Kanal 2 aufgezeichneten Signale werden dann um die jeweilige Differenz zwischen der Generatorsignalvorgabe bezüglich des Signalmaximums nmax,gen und der auf Kanal 1 für jede Einzelmessung bestimmte Position des Signalmaximums nmax,Kanal1 korrigiert. Das Ergebnis dieser Kompensation ist in Abb. 5.9(b) dargestellt und in Tabelle 5.1 zusammengefasst. Der Jitter-Effekt konnte durch die dargestellte Maßnahme deutlich reduziert werden. Die verbleibende Differenz zwischen dem Mittelwert der kompensierten Positionen des Signalmaximums und der Generatorvorgabe beträgt 0,0102 Abtastwerte, was einer akzeptablen zeitlichen Differenz von ca. 200 ps entspricht. Für die durchzuführenden Transmissionsmessungen sind 10 4 Mittelungen deutlich zu viel. Es werden nur N= 50 Mittelungen im Zeitbereich durchgeführt, um zeitlich kompakte Experimente zu erreichen. Durch diese Einschränkung erhöht sich die Unsicherheit der Jitter-Korrektur geringfügig. 80 5 Realisierung des inversen Messverfahrens Die Stichprobenstreuung für 50 Mittelungen liegt dann bei sσ nmax = 0,174 und die entsprechende Unsicherheit des Stichprobenmittelwertes bei unmax =snmax /√N=0,0246 [JCG08a]. Tabelle 5.1: Zusammenfassung der Merkmale des trigger jitters vor und nach der Kompensation. Arith. Mittelwert nmax Stichprobenstreuung sσ nmax Syst. Abweichung Vor Kompensation 103,1455 1,156 1,0776 Nach Kompensation 104,2681 0,1627 0,0102 Zählindex der Messung Abtastwert des Signalmaximums nmax 0 100 200 300 100 101 102 103 104 105 106 Kanal 1 Kanal 2 (a) Berechnete Position des Signalmaximums im Verlauf der ersten 300 Messungen. Komp. Abtastwert des Signalmaximums nmax,kom Rel. Häufigkeitsdichte 103.5 104 104.5 105 0 0.5 1 1.5 2 2.5 3 (b) Berechnete Position des Signalmaximums auf Kanal 2 nach Kompensation des trigger jitter-Merkmals von Kanal 1. Abbildung 5.9: Kompensation des trigger jitters auf Kanal 2 durch Ausnutzung des Jitter-Merkmals auf Kanal 1. 5.2.3 Messeffekt und Probekörperdimensionierung Der bei der Transmissionsmessung nach Abb. 5.2 genutzte Messeffekt wurde bereits von Rautenberg [Rau12] näher betrachtet und wird hier kurz als vereinfachtes Modell beschrieben, um anschließend eine Dimensionierung der Probekörper vorzunehmen. Dazu wird zunächst ein isotroper, elastischer, zylindrischer Stab angenommen, wie er bereits zur Materialparameterbestimmung nach Reynolds im Einsatz war [Rey53]. Für eine erste Analyse werden, trotz der eingeschränkten Geometrie, Freiraumwellen (Longitudinalwellen L mit Longitudinalwellengeschwindigkeit cL∈R und Transversalwellen T mit Transversalwellengeschwindigkeit cT∈R ) angesetzt. Von einem Longitudinalprüfkopf ausgehend, wird eine Longitudinalwelle L0 in den Probekörper eingekoppelt, siehe Abb. 5.10. An der Grenzfläche vom Probekörper zur Umgebungsluft kommt es nach dem Reflexionsgesetz der Akustik sin θ1 sin θ2=cL cT(5.7) 5.2 Konzeption und Funktionalität des Versuchsaufbaus 81 zur Reflexion einer Longitudinalwelle und einer Transversalwelle, wobei der Ausfallwinkel der Longitudinalwelle dem Einfallwinkel entspricht. Entsprechend ergibt sich, wie in Abb. 5.10 dargestellt, eine direkte Longitudinalwelle L0 0 durch den Prüfkörper und eine Transversalwelle T0 0 , welche unter dem Winkel θ2 durch den Prüfkörper läuft. An der gegenüberliegenden Grenzfläche zur Umgebungsluft wiederholt sich der Prozess der sogenannten Mode-Konversion von einer Transversalwelle nun in eine reflektierte Transversalwelle T0 1 und eine konvertierte Longitudinalwelle L0 1 . Dieser Vorgang kann sich dann im Probekörper in Abhängigkeit vom Winkel θ2 und der Länge l des Prüfkörpers n -Mal wiederholen. Am Ende des Prüfkörpers wird mit einem Longitudinalprüfkopf das Empfangssignal aufgezeichnet, welches sich zu einer Folge der direkten Longitudinalwelle L0 0 und der konvertierten Longitudinalwellen L0 nergibt. Der Eintreffzeitpunkt der direkten Longitudinalwelle L0 0 ergibt sich unter Berücksichtigung der Probenlänge lzu t0=l cL. (5.8) Die Zeitdifferenz ∆tzwischen dem Eintreffen der Wellen L0 0und L0 1kann bestimmt werden zu ∆t=Da cTcos θ2−Datan θ2 cL, (5.9) wobei Da∈R+ den (Außen-)Durchmesser der Probe beschreibt. Obige Gleichungen lassen sich nun umschreiben zu cL=l t0 cT=cL q1+cL∆t Da2, (5.10) welche genutzt werden können, um aus den gemessenen Empfangssignalen Startwerte, zur Lösung eines inversen Problems unter Berücksichtigung wellentheoretischer Aspekte im Vorwärtsmodell, zu schätzen. Das Strahlenmodell kann ebenso für eine hexagonale Materialsymmetrie angesetzt werden, siehe [Rau12], wobei eine Unterscheidung zwischen den Eigenschaften in der Vorzugsrichtung und der isotropen Ebene des Materials vorzunehmen ist. Allerdings stehen zur Lösung des entstehenden Gleichungssystems zu wenig Gleichungen zur Verfügung. Obige Ausführungen zum Messeffekt beziehen sich auf einen Strahlenmodellansatz, welcher aufgrund des Verhältnisses von Wellenlänge zum Probendurchmesser nur eingeschränkt Gültigkeit hat. Wird der Messeffekt unter wellentheoretischen Aspekten betrachtet, so zeigt sich, dass die beiden separierten Echos durch eine konstruktive Überlagerung von Signalanteilen der angeregten Moden entstehen. Im zeitlichen Bereich zwischen den Echos überlagern sie sich entsprechend destruktiv. Durch die räumliche Anregung auf der Stirnseite werden im Wellenleiter ausschließlich longitudinale Moden angeregt, davon allerdings entsprechend dem Frequenzbereich und den Materi- 82 5 Realisierung des inversen Messverfahrens Abbildung 5.10: Messeffekt dargestellt am zylindrischen Wellenleiter, Abb. nach [Rau12]. aleigenschaften der Probe sehr viele (größer 15), sodass möglichst viele Freiheitsgrade des Materials aktiviert werden und somit eine Sensitivität auf diese Parameter entsteht. Nähere Ausführungen zur Sensitivität werden in Kapitel 5.6 gegeben. Zur Probekörperdimensionierung wird im Folgenden kurz auf die verwendeten Kriterien zur Auslegung eingegangen, welche bereits von Rautenberg [Rau12] erarbeitet worden sind. Neben den oben benannten Signalbestandteilen, wird das experimentell bestimmte Empfangssignal ebenfalls Mehrfachechos enthalten, welche durch wiederholtes Laufen der Welle von einem zum anderen Ende des Prüfkörpers entstehen. Es wird durch die Dimensionierung der Probekörper erzielt, dass die Echos L0 0 und L0 1 möglichst gut separiert voneinander und ohne Überlagerung mit weiteren Echos vorliegen. Dazu lassen sich insgesamt drei Bedingungen formulieren: • Die Wanddicke Dw des Wellenleiters muss mindestens so groß sein, dass die entstehende Zeitdifferenz ∆t größer ist als die Pulsdauer Tp der Echos. Mit der Poisson-Zahl des isotropen Materials νiso ergibt sich für die minimale Wandstärke die Gleichung ∆t≥Tp⇒Dw,min =TpcLp1−2νiso . (5.11) • Die konvertierte Transversalwelle T0 0 muss die gegenüberliegende Wand erreichen können, um sich in die Longitudinalwelle L0 1 zu wandeln. Es ergibt sich für die Länge des Wellenleiters l>Dwtan θ2⇒lmin,1 =Dwp1−2νiso . (5.12) 5.3 Schallsender und -empfänger 83 • Das Mehrfachecho von L0 0 , welches nach dreimaligem Durchlaufen der Probe wieder am Empfänger sichtbar wird, darf sich nicht mit dem Echo L0 1 überlagern. Entsprechend ergibt sich für die Länge leine zweite Bedingung 2t0≤∆t+Tp⇒lmin,2 =1 2 Dw √1−2νiso +TpcL!, (5.13) wobei entsprechend die Bedingung mit der größeren Minimallänge anzuwenden ist. Neben den theoretischen Betrachtungen sind noch Aspekte wie die Verfügbarkeit der Probekörper, die Größe der aktiven Schallwandleroberfläche sowie die stark dämpfenden Eigenschaften der Proben zu berücksichtigen. Alle drei Aspekte zielen auf möglichst kleine Proben. Für die Auslegung der Probengeometrie wird im Folgenden eine Pulsdauer von Tp= 3 µ s zugrundegelegt. Zunächst wird die Wanddicke Dw ausgelegt. Dazu wird Gl. (5.11) für einen großen Wertebereich bezüglich der Poisson-Zahl und der Longitudinalwellengeschwindigkeit ausgewertet, siehe Abb. 5.11(a). Auf Basis einer Reihe von Literaturdaten bezüglich der Longitudinalund Transversalwellengeschwindigkeit von thermoplastischen Polymeren, kann der betrachtete Bereich zunächst eingegrenzt werden auf cL= [ 1250 m s−1 .. 3000 m s−1] und νiso = [ 0,25 .. 0,5 [ . In den Abbildungen 5.11(a) und 5.11(b) sind die genutzten Daten illustriert mit folgender Darstellungsweise: PMMA + , PP o, PTFE ∗ , LDPE × , HDPE 2 , PA ♦ , PEEK  , PMP F , PS  , POM  . In Anhang B.4 sind die Daten tabellarisch mit den entsprechenden Quellenangaben zusammengefasst. In Abb. 5.11(a) zeigt sich für den Punkt Dw,min( 3000 m s−1 ; 0,25 )≈ 6,36 mm die für diesen betrachteten Bereich zulässige kleinste Wandstärke. Im Folgenden soll daher die Wandstärke der verwendeten Proben etwa 6,5 mm betragen. Des Weiteren kann nun die Länge der Probe ausgelegt werden, siehe Abb. 5.11(b). Für Poisson-Zahlen nahe 0,5 steigt die minimal notwendige Länge der Probe stark an, sodass der Bereich bezüglich der Poisson-Zahl auf νiso = [ 0,25 .. 0,465 [ reduziert wird. Diese Einschränkung umfasst jedoch weiterhin alle hier aufgeführten Daten aus der Literatur. Für den Punkt lmin( 3000 m s−1 ; 0,465 )≈ 17,7 mm ergibt sich die minimal zulässige Probenlänge, sodass im Folgenden eine Länge von etwa 18 mm verwendet werden soll. 5.3 Schallsender und -empfänger 5.3.1 Anforderungen und gegenständlicher Entwurf Für das inverse Verfahren ist die Kenntnis über die Dynamik der Schallwandler (sowie der Verstärker, vlg. Kapitel 5.4) von zentraler Bedeutung. Das für den inversen Ansatz genutzte Empfangssignal beinhaltet die Eigenschaften der Sendeund Empfangselektronik, der Schallwandler und der zu charakterisierenden Probe. Aus diesem Grund muss das eingesetzte Messsystem hinreichend gut identifiziert sein, bevor belastbare Aussagen über die Eigenschaften von Probematerialien 84 5 Realisierung des inversen Messverfahrens Schallgeschw. cL/m s−1 Poisson-Zahl νiso 1500 2000 2500 3000 3500 4000 0 0.1 0.2 0.3 0.4 2 4 6 8 10 Dw/mm (a) Auslegung der Wandstärke des Probekörpers. Schallgeschw. cL/m s−1 Poisson-Zahl νiso 1500 2000 2500 3000 3500 4000 0 0.1 0.2 0.3 0.4 10 15 20 25 30 l/mm (b) Auslegung der Länge des Probekörpers. Abbildung 5.11: Auslegung der Probengeometrie. getätigt werden können. Hinsichtlich der Schallwandler ist die Rückwirkung der akustischen Schallkennimpedanz des Zielmediums auf deren dynamisches Verhalten zu berücksichtigen. Es ist also an dieser Stelle nicht zielführend, an einer Referenzprobe das dynamische Verhalten des Messsystems zu bestimmen. Vielmehr muss ein Modell des Schallwandlers aufgestellt und identifiziert werden, welches die Schallkennimpedanz des Zielmediums als freie Variable enthält. Neben recht komplexen 2D oder 3D FEM Simulationen [HRU09; WLB08; MBA06] existieren auch analytische Beschreibungsmodelle (KLM-Modell, Mason-Modell, usw.), die sich auf eine räumliche Dimension beschränken und somit in ihrer Gültigkeit eingeschränkt sind [LSW09; Mas48; KLM70]. Unter der Voraussetzung eines rotationssymmetrischen, nahezu reinen Dickenschwingers liefern diese jedoch hinreichend gute Ergebnisse. Gegenständlich werden die eingesetzten Schallwandler so ausgelegt, dass sie drei zentrale Anforderungen erfüllen: •Große Bandbreite der Wandler bei dennoch akzeptabler Signalamplitude • Gleichförmige Auslenkung bzw. geringe tangentiale mechanische Spannungskomponenten an der aktiven Schallwandlerfläche • Modellierung des Schallwandlers soll in einer räumlichen Dimension als Mason-Modell hinreichend gut möglich sein Diese Anforderungen werden aus konstruktiven Gesichtspunkten erfüllt, indem der Schallwandler folgende Eigenschaften aufweist: Als aktives piezoelektrisches Element wird ein 1-3 Piezokomposit verwendet. Dessen akustische Impedanz ( < 15 MRayl ) liegt deutlich unterhalb der Impedanz von gängigen Piezokeramiken (ca. 20 − 25 MRayl ), sodass keine akustischen Anpassschichten verwendet werden müssen, um den Wandler an die akustische Impedanz der Polymerproben anzupassen. 5.3 Schallsender und -empfänger 85 Elektrode PZT Epoxidharz Abbildung 5.12: Skizze eines 1-3 Piezokomposit in der dice and fill Ausführung. Abb. 5.12 zeigt ein Prinzipbild der verwendeten 1-3 Piezokomposite in der dice and fill Technik. Der Abstand der Bleizirkonat-Titanat-Stäbchen (PZT-Stäbchen) in einer Epoxidharz-Matrix ist so ausgelegt, dass der pitch von einem zum nächsten Stäbchen der halben Wellenlänge der auf der Platte ausbreitungsfähigen Grundmode der Plattenwelle entspricht. Die Wellenlänge wird an der Dickenresonanzfrequenz, also der späteren Betriebsfrequenz, ausgelegt. Dadurch überlagern sich die lateral auf der Komposit-Platte laufenden Wellen, welche durch die Lateralbewegungen der einzelnen Stäbchen erzeugt werden, destruktiv. Hingegen wirken in der Dickenbewegung alle Stäbchen konstruktiv zueinander. Dadurch ist es mit den 1-3 Piezokompositen möglich, nahezu ideale Dickenschwinger zu realisieren. Entsprechend sind die tangentialen Verschiebungen (radiale Richtung) eines aus diesem Material gefertigten Schallwandlers gering im Vergleich zu den Normal-Verschiebungen. Einhergehend kann aufgrund dieser Eigenschaft der (nahezu) reinen Dickenschwingung das aktive Element in nur einer Dimension modelliert werden. Um die Schwingeigenschaften der Piezokomposit-Scheibe möglichst nicht zu stören, wird auf die Nutzung eines Umkontakts seitens der Elektrode verzichtet. Die elektrische Kontaktierung der Piezokomposit-Scheibe erfolgt rückseitig durch einen dünnen flexiblen elektrischen Leiter. Vorderseitig erfolgt die elektrische Kontaktierung über eine Edelstahl-Schale mit einer Wanddicke von 75 µ m, vergleiche Abb. 5.13. Die Schale übernimmt dabei primär die Aufgabe einer Schutzschicht gegenüber mechanischen und chemischen äußeren Einflüssen und stellt sekundär die elektrische Verbindung zur vorderseitigen Elektrode (Masse-Potential) zur Verfügung. Das Masse-Potential wird dann von der Edelstahl-Schale über die Edelstahl-Hülse auf die Rückseite des Wandlers auf eine SMB-Steckverbindung weitergeführt. Als kritisch hat sich bei der Fertigung der Schallwandler der Klebstoff herausgestellt. Leitfähige Klebstoffe sind zumeist hochviskos und erreichen ihre Leitfähigkeit über kleine Silberkugeln. Dadurch ist die minimal erreichbare Klebschichtdicke recht groß ( > 40 µ m) [BRH11], was sich erschwerend auf eine planparallele Klebung auswirkt. Des Weiteren rückt die Dickenresonanzfrequenz des Schallwandlers mit größer werdender Schichtdicke 92 5 Realisierung des inversen Messverfahrens wobei Zi j die Komponenten der Impedanzmatrix aus Gleichung (5.26) bezeichnen. Einsetzen dieses Ausdrucks in die zweite und dritte Zeile von Gleichung (5.26) liefert ˜ F2=Z11 −Z2 12 Z11 +Zak,BA | {z } ˜ Z11 ˜ V2+ Z13 −Z12Z13 Z11 +Zak,BA! | {z } ˜ Z12 ˜ I(5.32) ˜ U= Z13 −Z12Z13 Z11 +Zak,BA! | {z } ˜ Z21 ˜ V2+Z33 −Z2 13 Z11 +Zak,BA | {z } ˜ Z12 ˜ I. (5.33) Die beiden Gleichungen können nun für den Sendeund für den Empfangsfall in Kettenmatrix-Form überführt werden. Für den Sendefall zeigt sich  ˜ F2 ˜ V2=˜ Z−1 12  ˜ Z11 ˜ Z11 ˜ Z22 −˜ Z2 12 1˜ Z22  | {z } AS  ˜ U −˜ I, (5.34) wohingegen sich für den Empfangsfall folgende Darstellung ergibt:  ˜ U ˜ I=˜ Z−1 12  ˜ Z22 ˜ Z11 ˜ Z22 −˜ Z2 12 1˜ Z11  | {z } AE=A−1 S  ˜ F2 −˜ V2. (5.35) In der folgenden gesamtheitlichen Betrachtung des Senderund Empfängermodells seien alle elektrischen Baugruppen m und alle akustischen Schichten n zusammengefasst zu jeweils einer Kettenmatrix Ael =Y m A(m) el (5.36) Aak =Y n A(n) ak . Wie in Abbildung 5.16 dargestellt, sind bei der senderseitigen Betrachtung die elektrischen und akustischen Kettenmatrizen zu invertieren. Des Weiteren ist zu beachten, dass im Sendefall das System mit der akustischen Last des Zielmediums abgeschlossen wird. Am elektrischen Eingangstor kann nun eine ideale Spannungsquelle angesetzt werden, vorausgesetzt, der Ausgangswiderstand der realen Quelle wird in Ael berücksichtigt. Außerdem wurden in Abbildung 5.16 zwei neue Blöcke AFT und ATF ergänzt, welche die Aufgabe der Transformation zwischen Kräften und Spannungen übernehmen mit AFT = A−10 0 1und ATF = A0 0 1. (5.37) 5.3 Schallsender und -empfänger 93 Durch Ausnutzung der Kettenmatrix-Formulierung kann nun direkt das Sendeübertragungsverhalten beschrieben werden durch  ˜σ0 ˜ V0=A−1 ak ·AFT ·AS·A−1 el | {z } MS · ˜ U0 −˜ I0. (5.38) Durch Einbeziehung der akustischen Impedanz des Medienabschlusses Zak,M mit ˜ V0=−˜σ0 Zak,M (5.39) ergibt sich der Sendefrequenzgang (in der Literatur oft Sendeübertragungsfaktor [LSW09]) MS∈C zu MS(iω)=˜σ0 ˜ U0 =Zak,M detMS MS,12 +MS,22Zak,M , (5.40) wobei MS,i j die Komponenten der Matrix MS beschreiben. Analog kann auch das Empfangsübertragungsverhalten berechnet werden. Zu beachten ist nun, dass auf der akustischen Seite das System angeregt wird und das elektrische Ausgangstor durch den elektrischen Eingangswiderstand des Signalverstärkers Zin,E ∈Cabgeschlossen wird. Es ergibt sich  ˜ U0 ˜ I0=Ael ·AE·ATF ·Aak | {z } ME · ˜σ0 −˜ V0mit ˜ I0=− ˜ U0 Zin,E (5.41) und erhält schließlich den Empfangsfrequenzgang (in der Literatur oft Empfangsübertragungsfaktor [LSW09]) ME∈Cmit ME(iω)=˜ U0 ˜σ0=Zin,E detME ME,12 +ME,22Zin,E . (5.42) Die Impulsantwort des modellierten Sendewandlers hS(t)∈R und Empfangswandlers hE(t)∈R berechnet sich nun durch inverse Fouriertransformation F−1{·} der jeweiligen Frequenzgänge. 5.3.3 Identifikation der Schallwandler Nachdem nun jeweils ein Modell für Sendeund Empfangswandler vorliegt, stellt sich die Frage nach der Parametrisierung der Modelle, um die real aufgebauten Schallwandler adäquat abzubilden. Die Identifikation der Modellparameter für jeden Schallwandler wird im Rahmen eines nichtlinearen inversen Problems realisiert. Als Messgröße dient die elektrische Eingangsimpedanz der Schallwandler Zin(iω)∈C, welche sich mit dem aufgestellten Modell bestimmen lässt zu Zin(iω) = ˜ U0 ˜ I0 =Zak,MMS,22 +MS,12 MS,11 +Zak,MMS,21 . (5.43) 94 5 Realisierung des inversen Messverfahrens Empfänger-Modell ˜ F˜ F1 ˜ V1 ˜ U0 ˜ U1 ˜ I1 ˜ I0 Zak,B Leitungsmodell und elektrisches Anpassnetzwerk Piezokomposit Backing Schutzschicht und Klebung ˜ F ˜ V ˜ F1 ˜ V1 ˜ U0 ˜ U1˜ I1 ˜ I0 Zak,B Zak,M Sender-Modell Leitungsmodell und elektrisches Anpassnetzwerk Piezokomposit Backing Schutzschicht und Klebung Kanal Eingangsimpedanz Empfangsverstärker ˜ V A−1 el A−1 ak Aea AFT ˜ V ˜σ ˜ V0 ˜σ0 −˜ V0 ˜σ0Aak ATF ˜ V ˜σAea Ael Zin,E Abbildung 5.16: Senderund Empfängermodell mittels elektrischer und akustischer 2-Tore sowie eines elektroakustischen 3-Tors. 5.3 Schallsender und -empfänger 95 Bei bekanntem Medienabschluss (z.B. Luft mit Zak,M ≈ 430 Rayl ) lässt sich somit jedes Element der Matrix MS in Zusammenhang mit der elektrischen Eingangsimpedanz bringen. Aufgrund der Reziprozität der Schallwandler [LSW09] (vgl. auch AE=A−1 S ) lässt sich aus den identifizierten Eigenschaften auf Basis der Sendeeigenschaft ebenfalls die Empfangsmatrix MEberechnen. Grundlegend ist die Parameteridentifikation in zwei Schritte unterteilt. In einem ersten Schritt werden die Eigenschaften des Piezokomposits bestimmt. Im zweiten Schritt werden dann die akustischen Schichten sowie das Backing ergänzt und somit der Schallwandler als Ganzes bei 21 ◦C identifiziert, wobei die Umgebungstemperatur durch Nutzung der in Kapitel 5.2.1 beschriebenen Thermokammer definiert eingestellt werden kann. Identifikation des Piezokomposits Zunächst wird das Modell des piezoelektrischen 1-3 Komposits identifiziert. Nach Gl. (5.26) sind fünf Modellparameter zu bestimmen: Zm , h33 , CS 0 , γ sowie tpc . Letzterer kann über eine direkte Messung der Scheibendicke mit Hilfe einer Messschraube bestimmt werden. Für die anderen Parameter sind folgende Zusammenhänge bekannt [LSW09; Rau12; ANS88; DIN76]: A=πD2 pc/4 (5.44) Zm=ρ˜cD TA(5.45) h33 =cD0 Tv tk2 Tρ εS 33 (5.46) CS 0=tpc AβS 33 (5.47) γ=ik=ω ˜cD T . (5.48) Gl. (5.44) beschreibt die Fläche der piezoelektrischen Scheibe in Abhängigkeit ihres Durchmessers Dpc , welcher direkt mit einer Messschraube gemessen wird. Zur Berechnung der weiteren Gleichungen seien die auftretenden Parameter im Folgenden benannt: •Dichte des Komposit-Materials: ρ∈R •Komplexe Schallgeschwindigkeit in Dickenrichtung: ˜cD T∈C •Realteil der komplexen Schallgeschwindigkeit: cD0 T∈R •Kopplungsfaktor der Dickenschwingung: kT∈R •Geklemmte Permittivität: εS 33 =βS 33−1 96 5 Realisierung des inversen Messverfahrens Startwertschätzung: Eine direkte Bestimmung dieser Parameter ist nicht möglich, jedoch können sie durch einige Näherungsgleichungen geschätzt werden, sodass Startparameter für eine nachfolgende Optimierung vorhanden sind. Als direkte Messgrößen sind zugänglich •Messung der Geometrie mit einer Messschraube ◦Dicke der Piezokomposit-Scheibe tpc ◦Durchmesser der Piezokomposit-Scheibe Dpc • Messung der elektrischen Eigenschaft mit einer Messbrücke (AIM-Cambridge LCR Databridge 401) ◦Freie Kapazität CT 0(bestimmt bei f=1 kHz) •Gravimetrische Messung mit einer Feinwaage (Sartorius LA 310 S) ◦Masse der Komposit-Scheibe inkl. Elektroden m • Analyse der elektrischen Eingangsimpedanz der Schallwandler Zin,mes( i ω)∈C mit dem Anritsu MS4630 Networkanalyzer ◦Minimalimpedanzfrequenz der Dickenresonanz fmin,d ◦Betragsimpedanz der Dickenresonanz Zmin,d =|Zin,mes(ω=2πfmin,d)| ◦Maximalimpedanzfrequenz der Dickenresonanz fmax,d ◦Minimalimpedanzfrequenz der Radialresonanz fmin,r ◦Betragsimpedanz der Radialresonanz Zmin,r =|Zin,mes(ω=2πfmin,r)| ◦Maximalimpedanzfrequenz der Radialresonanz fmax,r ◦ Minimalimpedanzfrequenz der ersten höheren Harmonischen der Radialresonanz fmin,r1 Die für jedes Piezokomposit bestimmten direkten Messgrößen sind im Anhang B.5 tabellarisch aufgeführt. Ebenfalls wird dort eine eindeutige Zuordnung der verschiedenen Schallwandler zu einer Identifikationsbezeichnung (Kennung) gegeben. Zunächst wird auf die Bestimmung der Dichte des Komposit-Materials ρ aus der Massenmessung mit ρeff= 4 m(πtpcD2 pc)−1 eingegangen. Um den Massenanteil der Kupfer-Zinn-Elektroden aus der Messung zu separieren, muss zunächst die Dicke der Elektrodenschicht te abgeschätzt werden. In der IEEE Norm 176 [ANS88] wird dazu eine Näherungsgleichung angegeben: te≈tpcρeff 2ρe ∆f/fmax,d −4k2 T/π2 1+∆f/fmax,d −4k2 T/π2, (5.49) 5.3 Schallsender und -empfänger 97 wobei ρe= 8500 kg m−3 die Dichte des Elektrodenmaterials und ∆f die Differenz zwischen Serienund Parallelresonanzfrequenz darstellt mit [GSC85] ∆f≈fmax,d −fmin,d q1+16 πfmin,dCT 0Zmin,d2. (5.50) Eine Schätzung für den Kopplungsfaktor kTist gegeben durch [ANS88] k2 T≈πfmin,d 2fmax,d tan π∆f 2fmax,d !. (5.51) Entsprechend kann die Dichte des Piezokomposit-Materials nun bestimmt werden zu ρ=tpcρeff−2teρe tpc −2te. (5.52) Für die Bestimmung der komplex-wertigen Schallgeschwindigkeit ˜cD T wird zunächst die komplexe Steifigkeit ˜ CD 33 betrachtet, welche in Speichermodul CD0 33 und Verlustmodul CD00 33 bzw. Speichermodul und mechanische Güte Qm∈Rgeschrieben werden kann ˜ CD 33 =CD0 33 +iCD00 33 =CD0 33 1+iQ−1 m. (5.53) Entsprechend kann die Schallgeschwindigkeit ˜cD T ausgedrückt werden über den Speichermodul und die mechanische Güte zu ˜cD T=sCD0 33 ρs1+i Qm=cD0 Ts1+i Qm. (5.54) Es wird an dieser Stelle also keine Frequenzabhängigkeit der Ausbreitungsgeschwindigkeit und Absorption betrachtet. Der Realteil der Schallgeschwindigkeit cD0 T kann über die Resonanzbedingung in Dickenrichtung abgeschätzt werden zu [BRH13] cD0 T=2tefffmax,d mit teff=tpc +2te(cD0 Tc−1 e−1). (5.55) In der Gleichung ist teff eine effektive Dicke, welche die Inhomogenität der Schallgeschwindigkeit in Dickenrichtung durch den Übergang vom Komposit-Material zum Elektrodenmaterial mit Schallgeschwindigkeit ce=5227 m/s berücksichtigt. Durch Umformulierung erhält man cD0 T=2fmax,d(tpc −2te) 1−4tefmax,dc−1 e . (5.56) 98 5 Realisierung des inversen Messverfahrens Um die mechanische Güte zu berechnen, braucht es noch weitere Zwischenwerte. Zunächst wird die freie Permittivität εT 33 bestimmt zu εT 33 =CT 0(tpc −2te) A, (5.57) mit deren Hilfe ein Schätzwert für die geklemmte Permittivität berechnen werden kann [ANS88; DIN76] εS 33 ≈(1−k2 T)(1−k2 P)εT 33. (5.58) Die Größe kP∈R beschreibt den planaren Kopplungsfaktor der Radialschwingung, welcher über die Poisson-Zahl des (als isotrop angenommenen) Piezokomposits νiso ∈R und die Radialresonanzeigenschaften geschätzt werden kann [ANS88] k2 P≈ ˜ J1(υ1) + νiso −1 ˜ J1(υ1)−2mit ˜ J1(υ1) = υ1J0(υ1) J1(υ1), (5.59) wobei Jn(x)die Besselfunktion n-ter Ordnung beschreibt. Das Argument der Besselfunktionen υ1 wird bestimmt zu υ1=η1 1+∆fr fmin,r !mit ∆fr≈fmax,r −fmin,r q1+16 πfmin,rCT 0Zmin,r2. (5.60) Der Parameter η1kann als Polynom bestimmt werden zu η1=0,4138 fmin,r1 fmin,r !2 −3,067 fmin,r1 fmin,r !+7,554 , (5.61) wobei das Polynom aus den diskreten Werten der Tabelle 12 in [ANS88] berechnet wurde. Schlussendlich können nun die geklemmte Kapazität CS 0 und die mechanische Güte Qm bestimmt werden zu CS 0=εS 33πD2 pc 4(tpc −2te)und Qm≈2πfmin,dCS 0Zmin,d ∆f2 f2 max,d  −1 . (5.62) Die so gewonnenen Startwerte für alle untersuchten Piezokomposite sind in Tabelle 5.2 zusammengefasst. Exemplarisch ist in Abb. 5.17 der simulierte und der gemessene Verlauf der Betragsimpedanz eines 1-3 Piezokomposits gezeigt, wobei die simulierten Werte unter Verwendung der geschätzten Startwerte p(0) berechnet wurden. Es zeigt sich bereits ein recht gutes Ergebnis, kleine Abweichungen sind gerade außerhalb des Dickenresonanzbereichs sichtbar. 5.3 Schallsender und -empfänger 99 Frequenz f/MHz Betragsimpedanz |Zin(iω)|/Ω 01234 100 102 Abbildung 5.17: Vergleich zwischen gemessener Betragsimpedanz und simulierter Betragsimpedanz unter Verwendung der Startwerte p(0) ; exemplarisch für den Wandler 1.5M −T (rot: Simulation, blau: Messung). Tabelle 5.2: Zusammenfassung der geschätzten Startwerte der verschiedenen Piezokomposite. Kennung ρ/kg m−3cD0 T/m s−1QmkTεS 33/ε0 750k-R 3696,6 3557,5 59,45 0,596 1040,9 750k-T 3672,7 3565,2 59,49 0,612 1000,8 1.0M-R 4042,5 3585,2 52,87 0,62 1217,7 1.0M-T 4038,2 3585,6 63,01 0,62 1205,2 1.5M-R 4086,9 3614,6 105 0,60 1078,2 1.5M-T 4103,2 3630,4 99,6 0,61 1066,2 2.0M-R 4187,2 3603,7 57,8 0,61 1073,2 2.0M-T 4172,9 3628,7 56,99 0,61 1076,6 2.5M-R 3895,2 3583,4 65,82 0,59 913,2 2.5M-T 3898,2 3569,6 50,2 0,58 947,5 100 5 Realisierung des inversen Messverfahrens Sensitivitätsanalyse: Sind Startwerte für das Modell des 1-3 Piezokomposits gegeben, können durch Formulierung des Parameterschätzproblems als nichtlineares inverses Problem, vgl. Kapitel 5.1.1, die gesuchten Parameter (Ursachen) im Sinne eines LS-Schätzers bestimmt werden. In diesem Abschnitt wird die Sensitivität der Modellfunktion auf die gesuchten Parameter betrachtet. Nur wenn sich durch Variation eines Modelleingangsparameters eine Änderung am Modellausgang feststellen lässt, ist eine Schätzung dieses Modelleingangsparameters möglich (Beobachtbarkeit). Bezüglich der zu minimierenden Kostenfunktion nach Gl. (5.2) ist es ausreichend, die Sensitivität auf Basis des Vorwärtsmodells zu überprüfen. Die partielle Ableitung der Kostenfunktion nach dem i-ten Element des Vektors pliefert ∂ ∂pi ε(p,y) = 2(f(p)−y)T∂f(p) ∂pi . (5.63) Im hier vorliegenden Fall ist als beobachtbare Wirkung die gemessene elektrische Betragsimpedanz |Zin,mes( i ω)| des 1-3 Piezokomposits gegeben. Die Messgröße liegt diskret an den äquidistanten Stützstellen n·∆ωvor, sodass die Beobachtungen in Vektorschreibweise als Zin,mes = [|Zin,mes(1)|, ...|Zin,mes(n)|, ...|Zin,mes(M)|](5.64) notiert werden. Entsprechend wird das Vorwärtsmodell als Zin,sim(p)∈RM bezeichnet. Der Ursachen-Vektor setzt sich zusammen zu p=hρ,cD0 T,kT,Qm,εS 33iT und beschreibt die Gesamtheit der Optimierungsvariablen. Gesucht ist also zunächst die Ableitung ∂Zin,sim(p)/∂pi∈RM . Zur Berechnung dieser Ableitung wird im Folgenden das Software-Paket ADiMat [BBL02] verwendet, welches die Möglichkeit des algorithmischen Differenzierens (AD) in MATLAB einbindet. Das AD ist in der Lage, bei vertretbarem Zeitbzw. Rechenaufwand numerisch exakte Ableitungsinformationen im Rahmen der Maschinengenauigkeit von bestehenden Computerprogrammen zu liefern [GW08]. ADiMat nutzt dazu den Ansatz der source code transformation, d.h. die MATLAB-Funktion wird an einen source code transformation server gesendet, dort umgeschrieben und als neues m-File zurückgesendet, welches bei Ausführung neben dem eigentlichen Funktionswert auch die Ableitungsinformation berechnet. Einige weitere Punkte zum AD sind in Anhang A.3 zusammengefasst. Im konkreten Fall, zur Berechnung von Ableitungen der Betragsimpedanz, ist eine Einschränkung zu beachten. Die komplexe Betragsfunktion ist formal nicht komplex differenzierbar, da sie nicht die Cauchy-Riemann-Differentialgleichungen erfüllt. Für den Fall, dass z∈R gilt und die komplexe Betragsfunktion f(z) = |z| als Abbildung f:R2→R2 interpretiert wird, gelingt es dennoch, die Ableitung zu bilden mit [PBC95] ∂Zin,sim(p) ∂pi =∂|Z(p)| ∂pi = Re Z(p)Re (∂Z(p) ∂pi)+Im Z(p)Im (∂Z(p) ∂pi) |Z(p)|. (5.65) 5.3 Schallsender und -empfänger 101 In der Quellcode-Umsetzung bedeutet dies, dass ADiMat bis zur Berechnung der komplex-wertigen Impedanz genutzt wird und dann die Ergebnisse nach Gl. (5.65) weiterverarbeitet werden. Bei der Quantifizierung der Sensitivitäten ist aufgrund der unterschiedlichen Größenordnungen der Optimierungsvariablen eine Normierung sinnvoll. Nach [Hil98] wird im Folgenden die skalierte Sensitivität dss, im Englischen dimensionless scaled sensitivity, definiert zu yjdsspi= ∂yj ∂pip∗p∗ i. (5.66) Sie beschreibt die skalierte Sensitivität der Ursache pi auf die Beobachtung yj . Die partielle Ableitung wird dabei an einem Arbeitspunkt p∗ ausgewertet und mit dem Wert der Ursache im Arbeitspunkt p∗ i gewichtet. Als Arbeitspunkt dient im Folgenden der berechnete StartwerteVektor p(0) . Abb. 5.18 zeigt, exemplarisch für den Wandler 1.5M-T, die so berechneten skalierten Sensitivitäten. Es zeigen sich folgende Eigenschaften des Vorwärtsmodells: • Die geklemmte Permittivität εS 33 zeigt eine konstante Sensitivität über den gesamten betrachteten Frequenzbereich. • Der Kopplungsfaktor kT ist nur lokal im Bereich der Resonanz (nicht der Antiresonanz) sensitiv. • Die Güte Qm ist sehr lokal im Betragsimpedanzmaximum und Betragsimpedanzminimum sensitiv. • Die Schallgeschwindigkeit cD0 T und die Dichte ρ sind im gesamten Resonanzbereich sensitiv. • Die Sensitivitäten der geklemmten Permittivität und der Güte sind geringer als die Sensitivitäten der übrigen Parameter. • Eine Änderung des Kopplungsfaktors hat lokal im Bereich der Antiresonanz den gleichen Verlauf der Sensitivität wie eine Änderung der Schallgeschwindigkeit cD0 T. Lösen des inversen Problems: Für die Wahl eines geeigneten Optimierungsalgorithmus ist eine Klassifizierung des vorliegenden Problems notwendig. Es handelt sich um ein nichtlineares Problem mit einer einzelnen Zielfunktion und mehreren Optimierungsvariablen. Die Optimierungsvariablen sind kontinuierlich mit Beschränkungen (box constraints). Die MATLAB-Routine fmincon bietet die Voraussetzungen, diesen Problem-Typ zu lösen. In näherer Ausführung wird ein interior point Algorithmus genutzt, welcher von Mathworks als besonders robust beschrieben wird [Mat15]. Die Jacobi-Matrix, berechnet durch ADiMat bzw. Gl. (5.65), wird dem Algorithmus direkt zur Verfügung gestellt. Auf die exakte Berechnung der Hesse-Matrix wird hingegen verzichtet und 108 5 Realisierung des inversen Messverfahrens der Betriebsfrequenz des Wandlers in Resonanz mit dessen statischer Kapazität ist. Problematisch ist, dass das L-Matching-Verfahren nur für eine Frequenz ausgelegt werden kann, sodass es sinnvoll für Systeme hoher Güte wie Luftschallwandler eingesetzt wird. Im Kontext dieser Arbeit hat sich für den realisierten Leistungsverstärker, vgl. Kapitel 5.4, herausgestellt, dass große Induktivitäten im Anpassnetzwerk zu Instabilitäten führen. Die notwendigen Parallelinduktivitäten würden im Bereich 1 mH für die 2,5 MHz -Wandler bis hin zu 24 mH für die 750 kHz -Wandler liegen. Aus diesem Grund wird im Folgenden auf die Verwendung von Induktivitäten im Anpassnetzwerk verzichtet. Zur Erhöhung der Bandbreite der Schallwandler können ebenfalls Ohm’sche Widerstände in Serie zum Wandler, vgl. elektrische Last in Abb. 5.24, geschaltet werden. Der Widerstand bewirkt eine Dämpfung des schwingenden Systems, wodurch sich dessen Bandbreite erhöhen lässt, siehe Abb. 5.21. Jedoch geht über diesen Widerstand auch elektrische Energie bei der Anregung des Wandlers verloren, sodass weniger elektrische in akustische Energie umgesetzt werden kann. Simulativ lässt sich mit den aufgestellten Modellen der Schallwandler der Maximalwert des SendeAmplitudengangs in Abhängigkeit vom gewählten Vorwiderstand berechnen. Daraus kann wiederum, durch Bezug auf den Maximalwert des Sende-Amplitudengangs ohne Vorwiderstand, eine Dämpfung ermittelt werden. Durch Linearisierung des Zusammenhangs lässt sich für kleine Vorwiderstände eine Empfindlichkeit der Dämpfung bezüglich des Widerstandswerts approximieren. Die Ohm’schen Widerstände werden entsprechend so dimensioniert, dass ein guter Kompromiss zwischen Bandbreite und Signalamplitude entsteht. Die berechneten Empfindlichkeiten der Dämpfung auf einen Vorwiderstand und der gewählte Widerstandswert der Sendewandler sind in Tabelle 5.6 zusammengefasst. Frequenz f/kHz Sende-Amplitudengang /dB 500 1000 1500 -12 -10 -8 -6 -4 -2 05Ω 10 Ω 15 Ω 20 Ω 25 Ω 30 Ω 35 Ω RV (a) Sende-Amplitudengang. -3 dB Bandbreite /kHz Vorwiderstand RV/Ω 5 10 15 20 25 30 35 480 500 520 540 560 580 Dämpfung /dB 5 10 15 20 25 30 35 -2 -1.5 -1 -0.5 0 0.5 (b) Dämpfung bezüglich System ohne Vorwiderstand und -3 dB Bandbreite. Abbildung 5.21: Betrachtung des Einflusses eines Ohm’schen Vorwiderstands auf den SendeAmplitudengang, beispielhaft für den Schallwandler mit der Kennung 750k-T. 5.3 Schallsender und -empfänger 109 Tabelle 5.6: Zusammenfassung der gewählten Vorwiderstände für die jeweiligen Sendewandler. Kennung 750k-T 1.0M-T 1.5M-T 2.0M-T 2.5M-T Empfindlichkeit /dB/Ω-0,07 -0,16 -0,39 -1,26 -1,51 Vorwiderstand / Ω 15 11,2 5,6 3 3 Empfangswandler Eine Erhöhung der Bandbreite der Empfangswandler wird im Folgenden nur für die Wandler 750k-R und 1.0M-R durchgeführt. Ansatz ist die Verwendung eines T-Netzwerks, wie in Abb. 5.22 dargestellt. Die Serienund Parallelkapazität der Schaltung verändern in ihrem Zusammenwirken das kapazitive Verhalten des Schallwandlers und somit dessen elektrischen Impedanzverlauf im niedrigen Frequenzbereich (Frequenzen kleiner Resonanzfrequenz). Hingegen verändert die Induktivität der Schaltung den elektrischen Impedanzverlauf im hohen Frequenzbereich (Frequenzen größer Resonanzfrequenz). Die verwendeten Bauteil-Dimensionierungen für die beiden Schallwandler sind in Tabelle 5.7 zusammengefasst. Der Effekt der Anpassnetzwerke wird anhand der aufgestellten Schallwandler-Modelle in Abb. 5.23 dargestellt. Tabelle 5.7: Zusammenfassung der Bauteil-Dimensionierung der verwendeten Anpassnetzwerke für die Empfangswandler. Kennung Cs/µFCp/pF Ls/µHRa/Ω 750k-R 2,2 330 9,4 10 1.0M-R 2,2 1 2,95 10 AnpassnetzwerkVerstärker Schallwandler CsLs Ra Cp Abbildung 5.22: Anpassnetzwerk für Empfangswandler. 110 5 Realisierung des inversen Messverfahrens Frequenz f/kHz Amplitudengang /dB 0 500 1000 1500 -20 -15 -10 -5 0 (a) Schallwandler 750k-R. Frequenz f/kHz Amplitudengang /dB 0 500 1000 1500 2000 2500 -20 -15 -10 -5 0 (b) Schallwandler 1.0M-R. Abbildung 5.23: Betrachtung des Einflusses des T-Netzwerks auf den Empfangs-Amplitudengang; ohne Anpassnetzwerk (blau) und mit Anpassnetzwerk (rot). 5.4 Sendeund Empfangselektronik Für die Realisierung des Messplatzes wird avisiert nur einen Leistungsverstärker für den Sendefall und nur einen Empfangsverstärker für den Empfangsfall für alle fünf realisierten SchallwandlerPaare zu nutzen. Dies ist vor allem für den Leistungsverstärker aufgrund der stark unterschiedlichen Lastimpedanzen bzw. der hohen kapazitiven Lasten (einige Nano-Farad) und der erforderlichen hohen Bandbreite eine Herausforderung14. Der im Folgenden genutzte Empfangsverstärker wurde aus der vorangegangenen Dissertationsarbeit von Dr. Jens Rautenberg übernommen [Rau12] und kann bei einer Verstärkung von 12 dB und 20 dB im gewünschten Frequenzbereich eingesetzt werden. Es handelt sich bei dem Verstärker um einen invertierenden Verstärker auf Basis eines LT1227 stromrückgekoppelten Operationsverstärkers. Die Eingangsimpedanz der Schaltung beträgt Zin,E =Rin + 1 /( i ωCin) mit Rin = 50 Ωund Cin =220 µF. 14 Ich möchte mich an dieser Stelle bei Dr. Jens Rautenberg und Christian Hoof für die Unterstützung bei dem Entwurf und der Realisierung der Verstärkerschaltungen bedanken. 5.4 Sendeund Empfangselektronik 111 5.4.1 Entwurf und Realisierung des Leistungsverstärkers Die Entwurfsspezifikationen des Sendeverstärkers können wie folgt zusammengefasst werden: •Bandbreite: ca. 100 kHz bis 5 MHz, •Ausgangsspannung (Spitze-Spitze) Uss =50 V, • hoher Ausgangsstrom zum Treiben von kapazitiven Lasten, vgl. Butterworth-van-DykeModell (BvD-Modell) eines piezoelektrischen Schallwandlers; die statischen Kapazitäten der Wandler variieren von ca. 2 nF für den 750 kHz -Wandler bis 10,4 nF für den 2,5 MHz -Wandler. Um diese Ziele zu erreichen, wird ein ausgangsseitig differentiell arbeitender Leistungsverstärker aufgebaut. Zur Spannungsverstärkung werden eine invertierende und eine nicht-invertierende Operationsverstärkerschaltung realisiert, vergleiche Abb. 5.24. Der Verstärkungsfaktor der Schaltung wird über die Widerstände R3 , R4 , R5 und R6 bestimmt ( V= 8). Die Operationsverstärker vom Typ AD817 können jedoch nicht bei den geforderten Lasten arbeiten. Um den notwendigen Ausgangsstrom zu erreichen, werden parallel sogenannte Slave-Amplifier betrieben. Jedem Spannungsverstärker (Master) werden jeweils drei Slave-Amplifier zugeordnet, welche als Impedanzwandler (unity gain buffer) betrieben werden. Das Konzept zur Ausgangsstromerhöhung durch Parallelschaltung von Operationsverstärkern ist [Ap12; Ap13] entnommen. Als Slave-Amplifier werden Operationsverstärker vom Typ LT1210 verwendet, welche sich durch hohe Ausgangsströme auszeichnen ( Imax = 1,1 A). Zur Stabilisierung der Stromaufteilung an den Ausgängen der LT1210 werden 1 Ω Widerstände ( R1 ) verwendet. Der LT1210 ist mit einem Kompensationsnetzwerk bei kapazitiver Belastung zu betreiben. Durch Verwendung einer Kapazität mit CC= 1 nF zwischen dem COMP-Pin und dem Ausgang sowie einem Widerstand R2= 820 Ωim Rückkopplungszweig zwischen invertierendem Eingang und Ausgang werden Überhöhungen im Amplitudengang vermieden. Der Ausgang des Leistungsverstärkers darf in dieser Bauform nicht mit dem Massepotential verbunden werden. Dies ist für die Nutzung von Schallwandlern als Last kein Problem, solange die Schallwandler nicht mit Masse-führenden Komponenten elektrisch leitend in Kontakt treten. Bei der Verwendung elektrisch leitender Probekörper bestünde die Gefahr, dass ein leitfähiger Kontakt vom Sendewandler zum Empfangswandler und von diesem auf das Massepotential des Empfangsverstärkers entsteht. Um dieses Problem zu umgehen, wird ausgangsseitig ein Übertrager genutzt, um die Sekundärseite galvanisch von der Primärseite zu trennen. Problematisch bei der Suche nach einem geeigneten Übertrager sind die zu übertragende Leistung und der gewünschte Frequenzbereich. Die geforderten Eigenschaften werden durch den Übertrager Cooltronics VPH50155-R erfüllt. Primärsowie sekundärseitig können jeweils drei Induktivitäten parallel oder in Reihe geschaltet werden, wodurch sich auch das Übersetzungsverhältnis ergibt. Da der Übertrager ausschließlich zur galvanischen Trennung genutzt werden soll, werden die jeweils drei Primär- 112 5 Realisierung des inversen Messverfahrens + - + - + - + - + - + - + - + - R2 Uin R3 R4 C1 C2 CC 1k R1 R2 CC R1 R2 CC R1 1kR5 R6 R2 CC R1 R2 CC R1 R2 CC R1 Uout RV C0,BvD R1,BvD C1,BvD L1,BvD 1:1 Abbildung 5.24: Schaltbild Leistungsverstärker mit elektrischem Lastmodell (BvD). 5.4 Sendeund Empfangselektronik 113 und Sekundärinduktivitäten parallel geschaltet. Aus dem Datenblatt [Coo15] lassen sich folgende Kenngrößen des Übertragers ablesen: •Primärinduktivität pro Windung Lbase =22,3 µH •Sekundärinduktivität pro Windung Lbase =22,3 µH •Streuinduktivität pro Windung Lσ base =0,225 µH •Sättigungsstromstärke pro Windung Isat,base =1,05 A •Verlustwiderstand pro Windung Rbase =0,0711 Ω Durch die Parallelschaltung der Windungen ergeben sich die Primärund Sekundärinduktivität zu L1=L2=Lbase . Des Weiteren ergibt sich eine maximale Stromstärke Imax,sat , welche zu einer Sättigung des Kernmaterials von ca. 30% führt: Imax,sat =Isat,base32=28,35 A. (5.68) Es stellt sich die Frage nach der Bandbreite des Übertragers. Dazu müssen zunächst der Streufaktor σtrans und der Kopplungsfaktor ktrans bestimmt werden zu σtrans =Lσ baseL−1 base =0,0105 und ktrans =p1−σtrans =0,9947. (5.69) Die Übertragungsfunktion (in der Laplace-Transformierten) eines mit der elektrischen Last ZL(s) belasteten Übertragers kann beschrieben werden als [Dos14; Wei07] GU(s) = sMZL(s) s2(L1L2M2) + s(L2R1+L1R2+L1ZL(s)) + R1R2+R1ZL(s), (5.70) mit M=ktrans √L1L2 sowie R1=R2= 3 Rbase . Die elektrische Last kann wiederum beschrieben werden durch das BvD-Modell (Schallwandler) in Serie zu einem Vorwiderstand mit ZL(s) = RV+s2L1,BvDC1,BvD +sC1,BvDR1,BvD +1 s3L1,BvDC1,BvDC0,BvD +s2C1,BvDC0,BvDR1,BvD +s(C0,BvD) + C1,BvD . (5.71) Um die Grenzen des Übertragers abzuschätzen, reicht es aus, die Übertragungsfunktionen für die Last-Bedingungen mit dem Schallwandler 750k-T und dem Schallwandler 2.5M-T zu betrachten, siehe Abb. 5.25. Die Parameter der BvD-Modelle können invers bestimmt werden, siehe Anhang B.8. In Abb. 5.25(a) ist zu erkennen, dass der niedrige Frequenzbereich mit den hohen Impedanzen als Last keine Probleme verursacht. Hingegen ist in Abb. 5.25(b) zu erkennen, dass im Bereich der Betriebsfrequenz des Schallwandlers 2.5M-T der linear-phasige Bereich verlassen wurde und das gekoppelte System aus Schallwandler und Übertrager deutlich in der Überhöhung des Amplitudengangs betrieben wird. Außerdem ist hier eine klare Rückwirkung der Last zu erkennen. 114 5 Realisierung des inversen Messverfahrens Die Forderung nach 5 MHz Bandbreite kann demnach mit dem Übertrager nicht ganz erfüllt werden, dennoch stellt die Nutzung des Übertragers einen guten Kompromiss zwischen Sicherheit im Betrieb und Bandbreite dar. Frequenz f/Hz Amplitudengang /dB 104106108 -20 -10 0 10 Phasengang /◦ -100 -50 0 50 (a) Last: 750k-T. Frequenz f/Hz Amplitudengang /dB 104106108 -40 -20 0 20 Phasengang /◦ -200 -100 0 100 (b) Last: 2.5M-T. Abbildung 5.25: Bode-Diagramm des verwendeten Übertragers unter zwei verschiedenen Lastbedingungen. 5.4.2 Identifikation Die dynamischen Eigenschaften der beiden realisierten Verstärker sollen im Folgenden als IIR-Filter (infinite impulse response filter) modelliert und deren Filterkoeffizienten identifiziert werden. Zur Identifikation wird jeweils ein Testsignal auf den Verstärkereingang gegeben und das Antwortverhalten am Verstärkerausgang beurteilt. Als Testsignal wird ein Rauschsignal (random binary noise signal: RBS-Signal) verwendet, welches mit einer Grenzfrequenz von 10 MHz tiefpassgefiltert wurde (Tiefpass 2ter Ordnung vom Typ Butterworth). Zur Identifikation der Filterkoeffizienten wird der Steiglitz-McBride-Iterationsalgorithmus verwendet [SM65], wie er von MATLAB zur Verfügung gestellt wird. Zur Kontrolle der Schätzergebnisse sowie zur Abschätzung des Zählerund Nennerpolynomgrads wird zusätzlich eine nicht-parametrische Schätzung des Amplitudenund Phasengangs auf Basis der Kreuzund Autoleistungsdichtespektren berechnet. Ebenfalls wird neben den Testsignalen auch ein Validierungssignal aufgenommen, um die auf Basis der Testsignale berechneten Filterkoeffizienten zu überprüfen. Empfangsverstärker Der Empfangsverstärker kann unabhängig vom genutzten Schallwandler identifiziert werden, jedoch muss für jeden Verstärkungsfaktor, 12 dB und 20 dB , eine Identifikation der IIR-Filterkoeffizienten 5.4 Sendeund Empfangselektronik 115 durchgeführt werden. Eine ausführliche Betrachtung der einzelnen Stationen zur Identifikation und Validierung sei hier für den Fall der Verstärkung von 12 dB vorgenommen. Die Identifikation für 20 dB erfolgt analog. In Abb. 5.26(a) sind das Testsignal am Systemeingang und die Antwort des Systems im Zeitbereich dargestellt. Abb. 5.26(b) zeigt eine nicht-parametrische Schätzung des Amplitudengangs und des Phasengangs, wie sie sich aus der Schätzung mittels Kreuzleistungsdichte Sxy(iω)und Autoleistungsdichte Sxx(iω)ergeben: ˆ Gnonparam(iω) = Sxy(iω) Sxx(iω). (5.72) Es zeigt sich der Verlauf einer Übertragungsfunktion geringer Ordnung. Der Frequenzbereich größer 10 MHz ist aufgrund der geringen spektralen Anregung in diesem Bereich von einer kleinen Kohärenz geprägt und somit aus der Betrachtung auszuschließen. Zur parametrischen Systemidentifikation mit dem Steiglitz-McBride-Iterationsalgorithmus werden ein Nennerpolynomgrad von zwei und ein Zählerpolynomgrad von eins gewählt. Das Ergebnis dieser Identifikation ist in Abb. 5.26(c) dargestellt. Aus der parametrischen Schätzung lässt sich die − 3 dB -Grenzfrequenz zu 9,2 MHz bestimmen. Abb. 5.26(d) zeigt das verwendete Validierungssignal, eine Folge von Gaußpulsen steigender Mittenfrequenz von 750 kHz bis 2,5 MHz . Sowohl in Abb. 5.26(a) als auch in Abb. 5.26(d) zeigen die gemessenen Ausgangssignale (blau) und die auf Basis des identifizierten IIR-Filters berechneten Ausgangssignale (rot) eine gute Übereinstimmung. Zur quantitativen Bewertung der Qualität des identifizierten Filters wird auf Basis der Validierungssignale die Normierte RMS Abweichung (Normalized Root-Mean-Square Deviation: NRMSD) berechnet. Auf Basis der NRMSD kann iterativ eine Adaption der Vorgaben für die Polynomgrade des Zählers und des Nenners erfolgen. Als Definition des NRMSD wird NRMSD =100% · 1−kymess −ysimk2 kymessk2!(5.73) gewählt. Wird die Verstärkung des Empfangsverstärkers auf 20 dB erhöht, so verringert sich die − 3 dB -Grenzfrequenz auf ca. 5,3 MHz . Das Ergebnis der parametrischen Schätzung sowie das Einund Ausgangssignal zur Validierung der Schätzung sind in Abb. 5.27 dargestellt. Die identifizierten Filterkoeffizienten sind im Anhang B.9 aufgeführt. Leistungsverstärker Die dynamischen Eigenschaften des Leistungsverstärkers sind abhängig von der elektrischen Last (BNC-Leitung, Vorwiderstand und Schallwandler), d.h. die Identifikation des Leistungsverstärkers muss unter Berücksichtigung der Rückwirkung für jeden Sendewandler separat erfolgen. Als Evaluierungssignal wird weiterhin das gefilterte RBS-Signal verwendet. Zur Verifikation wird 116 5 Realisierung des inversen Messverfahrens Zeit t/µs Eingangssignal /V 0 2 4 6 8 10 12 -0.1 0 0.1 Zeit t/µs Ausgangssignal /V 0 2 4 6 8 10 12 -0.2 0 0.2 (a) Eingangsund Ausgangssignal zur Identifikation. Frequenz f/Hz Amplitudengang /dB 105106107108 0 5 10 15 20 Frequenz f/Hz Phasengang /° 105106107108 0 50 100 150 200 250 (b) Nicht-parametrische Schätzung. Frequenz f/Hz Amplitudengang /dB 104105106107108 0 5 10 15 Frequenz f/Hz Phasengang /° 104105106107108 0 100 200 (c) Parametrische Schätzung. Zeit t/µs Eingangssignal /V 0 5 10 15 20 -0.05 0 0.05 0.1 Zeit t/µs Ausgangssignal /V 0 5 10 15 20 -0.4 -0.2 0 0.2 (d) Eingangsund Ausgangssignal zur Validierung. NRMSD =96,04%. Abbildung 5.26: Systemidentifikation des Empfangsverstärkers bei V=12 dB Verstärkung. Frequenz f/Hz Amplitudengang /dB 104105106107108 -20 -10 0 10 20 30 Frequenz f/Hz Phasengang /° 104105106107108 0 50 100 150 200 (a) Parametrische Schätzung. Zeit t/µs Eingangssignal /V 0 5 10 15 20 -0.05 0 0.05 0.1 Zeit t/µs Ausgangssignal /V 0 5 10 15 20 -1 -0.5 0 0.5 (b) Eingangsund Ausgangssignal zur Validierung. NRMSD =96,33%. Abbildung 5.27: Systemidentifikation des Empfangsverstärkers bei V=20 dB Verstärkung. 5.4 Sendeund Empfangselektronik 117 jeweils ein Gaußpuls verwendet, dessen Mittenfrequenz zur Betriebsfrequenz des Sendewandlers passt. Die relative Bandbreite wird konstant zu 110% gewählt. Es stellt sich zunächst wieder die Frage nach dem zu verwendenden Zählerund Nennergrad zur Systemidentifikation. Wird die Verstärkerschaltung ohne Übertrager als Verzögerungsglied erster Ordnung approximiert und werden die Übertragungsfunktion des Übertragers aus Gl. (5.70) und die Last nach Gl. (5.71) einbezogen, so ergibt sich insgesamt ein Zähler der Ordnung sieben und ein Nenner der Ordnung neun. Tatsächlich lässt sich mit diesen Werten der Leistungsverstärker mit den Lasten 1.5M-T, 2.0M-T und 2.5M-T gut identifizieren (deutliche Rückwirkung auf den Verstärker). Für die Lasten 750k-T (Zählergrad 2, Nennergrad 3) und 1.0M-T (Zählergrad 4, Nennergrad 5) können die gewählten Zählerund Nennergrade reduziert werden, da hier die Rückwirkung der Last gering ausfällt. Abb. 5.28 bis Abb. 5.32 zeigen zusammenfassend die Ergebnisse der parametrischen Schätzung sowie das entsprechende gemessene (blau) und simulierte (rot) Validierungssignal. Die identifizierten Filterkoeffizienten sind im Anhang B.9 aufgeführt. Es zeigt sich, dass die Bandbreite des Leistungsverstärkers gerade ausreichend ist. Für den Schallwandler 2.5M-T mit einer Mittenfrequenz von 2,5 MHz liegt der Arbeitsbereich des Schallwandlers im Bereich der Grenzfrequenz und somit nicht mehr im linear-phasigen Arbeitsbereich, wodurch der Sendepuls deutlich verzerrt wird. Dieses Verhalten wurde in den theoretischen Überlegungen zum Übertrager bereits postuliert. Frequenz f/Hz Amplitudengang /dB 104105106107 0 10 20 Frequenz f/Hz Phasengang /◦ 104105106107 0 100 200 (a) Ergebnis der parametrischen Schätzung. Zeit t/µs Ausgangssignale /V 0123456 -30 -20 -10 0 10 (b) Validierungssignal: Messung (blau), Simulation (rot). NRMSD =98,82%. Abbildung 5.28: Systemidentifikation des Leistungsverstärkers bei Belastung mit Schallwandler 750k-T. 124 5 Realisierung des inversen Messverfahrens Anhand einer PEEK Probe mit 20 mm Länge wird die Reproduzierbarkeit der Koppelschichtdicken bewertet. Dazu werden mit dem Schallwandlerpaar 750k vier Transmissionsmessungen durch diese Probe durchgeführt, wobei nach jeder Einzelmessung das Koppelmittel vollständig zu entfernen und wieder neu aufzutragen ist. Ausgewertet wird die Laufzeit durch die Probe, wobei der Betrag der analytischen Korrelationsfunktion (siehe Anhang B.11) zwischen gemessenem Empfangssignal und simuliertem Empfangssignal ohne Proben-Einfluss15, siehe Abb. B.10, verwendet wird. Die Stichprobenstreuung der Laufzeitmessungen liegt dabei bei sσ=0,0074 µs. Durch die geringe Anzahl der Messungen ergibt sich der Erweiterungsfaktor zur Bewertung eines 95%-Vertrauensintervalls zu 3,18 und die erweiterte Unsicherheit zu U95% = 0,0236 µ s [JCG08a], was bei der verwendeten Abtastrate von 50 MHz etwa einem Abtastwert entspricht. (a) Siebdruck-Gitter. (b) Siebdruck. (c) Realisierung der Transmissionsstrecke. Abbildung 5.36: Realisierung des experimentellen Aufbaus zur Transmissionsmessung. Totzeitschätzung: Zur Schätzung der Totzeit (inkl. Laufzeit in den Koppelschichten) werden insgesamt sechs Stäbe aus PEEK mit konstantem Durchmesser von 20 mm aber variierender Länge li von 5 mm bis 40 mm genutzt. Für jede Länge wird eine Laufzeit Ti durch die Probe ermittelt und auf Basis der gemessenen Daten eine Kennlinienapproximation mit linearem Funktionsansatz im Sinne kleiner Fehlerquadrate ermittelt. Eine Extrapolation auf eine Probe mit virtueller Länge 0 m zeigt dann die gesuchte Totzeit Ttot des Messsystems an. 15 Im Folgenden sollen die mit den Impulsantworten von Sendeverstärker, Sendewandler, Empfangswandler und Empfangsverstärker gefalteten (und später auch um die hier zu ermittelnde Totzeit korrigierten) Sendesignale als Signalprototypen yp(t) bezeichnet werden. Sie bezeichnen die modellierte Signalformen, wie sie ohne Einfluss einer Probe im Signalpfad als Empfangssignal aufgenommen werden würden. 5.5 Analyse des Messsystems im Zeitbereich 125 Die Koeffizienten der Laufzeitfunktion T(l) = Ttot +mtotTergeben sich zu T=1 6X i Ti;l=1 6X i li T l =1 6X i Tili;l2=1 6X i l2 i Ttot =T−mtotl;mtot =l T l T l2l2. (5.75) Zur Betrachtung der Unsicherheit der Totzeitschätzung nach [JCG08a] wird die Unsicherheit der Längenmessung sowie der Laufzeitbestimmung im Sinne von Typ B-Unsicherheiten berücksichtigt. Der zur Längenbestimmung verwendeten Messschraube wird laut Datenblatt eine ’Genauigkeit’ von ± 3 µ m zugeordnet. Im Folgenden wird entsprechend die Varianz u2 µm,TypB = 3 2/ 3 (µ m )2 verwendet. Aus den Experimenten zur Darstellung der Reproduzierbarkeit ist eine Unsicherheit der Laufzeitbestimmung im Sinne von Typ B verwendbar mit u2 T,TypB = ( 1,2 · 0,0074 µ s )2 . Die kombinierte Standardunsicherheit uTtot c der Totzeitbestimmung ergibt sich nach Gl. (B.48). Der beste Schätzwert der Totzeitbestimmung wird mit Ttot gekennzeichnet und ergibt sich, da ausschließlich Typ B-Unsicherheiten betrachtet wurden, direkt durch Auswertung von Gl. (5.75). Die erweiterte Unsicherheit zur Berechnung des 95%-Vertrauensintervalls ergibt sich zu U95% Ttot = 1,96 uTtot c . Die ermittelten Werte sind Tabelle 5.8 zu entnehmen. Tabelle 5.8: Zusammenfassung der geschätzten Totzeiten sowie deren zugehörige erweiterte Unsicherheit für die fünf Wandler-Paare. Kennung 750k 1.0M 1.5M 2.0M 2.5M Ttot /µs 0,1981 0,1177 0,1423 0,0988 0,0536 uTtot c/µs 0,0071 0,0071 0,0071 0,0071 0,0071 U95% Ttot /µs 0,014 0,014 0,014 0,014 0,014 Die ermittelten Totzeiten beinhalten die Unsicherheit in den modellierten akustischen Systemen, die Totzeit der elektronischen Geräte, die Koppelschichten sowie die Unsicherheit der Korrelationsanalyse. Aus diesem Grund können die ermittelten Totzeiten von Wandler-Paar zu Wandler-Paar variieren. Die zugeordneten Unsicherheiten der Schätzwerte variieren nur im Nanosekundenbereich, was auf die implizite Linearisierung der Laufzeitfunktion im Punkt des besten Schätzwerts zur Berechnung der kombinierten Standardunsicherheit zurückzuführen ist. Der Punkt des besten Schätzwerts wird bezüglich der Längenangaben nicht verändert, da für jede Todeszeitbestimmung die gleichen Stäbe genutzt werden. Die mit den Schallwandler-Paaren gemessenen absoluten Laufzeiten unterscheiden sich hingegen geringfügig aufgrund der Frequenzabhängigkeit der Schallgeschwindigkeit. 126 5 Realisierung des inversen Messverfahrens Die ermittelten Totzeiten werden in der Vorwärtssimulation berücksichtigt. Die Verschiebung der Signale um fraktionale Teile der Abtastzeitpunkte wird durch ein Thiran-Filter realisiert, wie es im MATLAB implementiert ist (thiran) [LVK96]. Probenlänge l/mm Laufzeit T/µs 0 10 20 30 40 0 5 10 15 (a) Kennlinienapproximation. Probenlänge l/mm Laufzeit T/µs 0 0.1 0.2 0.3 0.4 0.5 0.15 0.2 0.25 0.3 0.35 0.4 (b) 95%-Unsicherheitsintervall der Approximation. Abbildung 5.37: Kennlinienapproximation der Laufzeit als Funktion der Probenlänge für Schallwandler-Paar 750k. Messdaten (rot), approximierte Kennlinie (blau). 5.6 Sensitivitätsanalyse Im Folgenden wird näher untersucht, inwiefern sich eine Änderung der gesuchten Materialparameter in einer Änderung des beobachtbaren Empfangssignals manifestiert. Wird ein viskoelastisches transversal isotropes Material mit dem fraktionalen Zener-Modell beschrieben, ergeben sich nach Gl. (4.55) vier linear unabhängige Eigenwerte und ein Eigen-Winkel zur Beschreibung des ’elastischen’ Parts. Hinzu kommen für jeden Eigenwert jeweils drei Parameter zur Beschreibung des ’viskoelastischen’ Verhaltens der entsprechenden Eigenbewegung. In Summe sind dies 17 unbekannte Parameter. In Kapitel 4.2.3 wurde bereits gezeigt, dass für die hier betrachteten Polymere der Fall der Quasi-Isotropie angenommen werden kann, sodass der Eigen-Winkel den Wert ωe= 125,27 ◦ annimmt. Des Weiteren wurde eine empirische Näherung für µL diskutiert, siehe Gl. (4.47), sodass das Materialmodell um bis zu zwei ’elastische’ Parameter reduziert werden kann. Die verbleibende Anzahl der Unbekannten ist dennoch beachtlich hoch, sodass in diesem Abschnitt eine Analyse der Sensitivität der (simulierten) Empfangssignale auf die unbekannten Materialparameter durchgeführt wird. Insbesondere wird dabei die Sensitivität der eigenbewegungsbezogenen Verlustmechanismen untersucht. Neben der skalierten Sensitivität nach Gl. (5.66) wird in diesem Abschnitt eine weitere Größe eingeführt, welche den Informationsgehalt einer Menge von Beobachtungen y bezüglich einer 5.6 Sensitivitätsanalyse 127 Ursache pi beschreibt. Diese Größe kann nach Hill [Hil98] als RMS-Wert der skalierten Sensitivität beschrieben werden zu (composite scaled sensitivity) ycsspi=v u t1 N N X j=1 (yjdsspi)2. (5.76) Als Beobachtungen y gelten die diskreten Werte der simulierten Transmissionsmessung, als Ursachen werden die gesuchten Materialparameter betrachtet. Die skalierte Sensitivität kann in Hinblick auf die Simulation der Transmissionsmessung nicht analytisch und auch nicht über AD berechnet werden. Daher wird im Folgenden auf die Bildung von Ableitungsapproximationen mittels zentralem Differenzenquotienten zurückgegriffen. Die Schrittweite zur Bildung der finiten Differenz wird zu ±1% der jeweiligen betrachteten Ursache angenommen. 5.6.1 Betrachtungen am elastischen transversal isotropen Material Als exemplarischer Arbeitspunkt gelte im Folgenden ein elastisches isotropes Material mit den Eigenschaften cL= 2300 m/s , cT= 1100 m/s und ρ= 1150 kg m−3 . Durch elementare Umformungen entspricht dies einem Kompressionsmodul von K≈ 4,23 GPa und einem Schermodul von µ≈ 1,39 GPa bzw. einem Elastizitätsmodul E≈ 3,76 GPa und einer Poisson-Zahl von νiso = 0,3517 [Mec08]. Ausgehend von dem isotropen Material wird durch Variation der unabhängigen elastischen Materialparameter (Ursachen) p= [EL , ET , νL , νT , µL]T ein Material hexagonaler Symmetrie erzeugt. Die Geometrie des Wellenleiters sei gegeben durch ri= 3 mm , Dw= 6,5 mm und l= 17,5 mm . Als zeitliche Anregung wird der Signalprototyp des 750k-Schallwandler-Paars genutzt. Da das elastische Material keiner Frequenzabhängigkeit unterliegt, steht diese spektral limitierte Betrachtung qualitativ für einen größeren spektralen Bereich. Abb. 5.38 zeigt die skalierten Sensitivitäten exemplarisch für EL und µL . Während EL über den gesamten Signalverlauf eine hohe Sensitivität zeigt, beschränkt sich der Einflussbereich von µL auf den Signalbereich, welcher nach dem Strahlenmodell von Reynolds mit den Anteilen der Transversalwelle korrespondiert. Werden die skalierten Sensitivitäten auch für die der anderen elastischen Materialparameter berechnet und zur composite scaled sensitivity umgerechnet, ergibt sich das in Abb. 5.39(a) dargestellte Verhältnis von hohem zu niedrigem Informationsgehalt der simulierten Zeitsignal-Beobachtungen. Die dargestellten Werte sind dabei auf den größten Informationsgehalt normiert. Wie in Abschnitt 4.2.3 diskutiert, kann für die betrachteten thermoplastischen Polymere der Fall der Quasi-Isotropie angenommen werden, sodass ein Materialparameter aus der Betrachtung entfällt, siehe Gl. (4.42). Es zeigt sich an dieser Stelle als sinnvoll, die hohe Sensitivität von νL auf die Materialparameter EL , ET und νT zu verteilen und somit im Folgenden νL als Funktion von EL , ET und νT auszudrücken. Der unter dieser Annahme berechnete Informationsgehalt ist in Abb. 5.39(b) dargestellt. Es zeigt sich eine leicht verbesserte Ausgewogenheit zwischen den 128 5 Realisierung des inversen Messverfahrens Informationsgehalten, mit Ausnahme des Informationsgehalts bezüglich µL . Dies ist nicht weiter verwunderlich, da sich µL auf die Scherung in der (θ , z) und (r , z) -Ebene bezieht und somit einen deutlichen Einfluss auf Torsionsmoden ausübt, welche zwar durch die Anisotropie nicht mehr entkoppelt von den Longitudinalmoden vorliegen, aber dennoch eine geringe Rolle in der Ausbreitung der Quasi-Longitudinalmoden einnehmen. Dieser Umstand zeigt, dass es bei der Lösung des inversen Problems bezüglich µL zielführender ist, diesen Materialparameter nach Gl. (4.47) zu schätzen, anstatt ihn als Optimierungsvariable zu betrachten. Zeit t/µs Empfangssignal 10 15 20 25 -0.5 0 0.5 1 (a) Variation von ELum ±1%. Zeit t/µs Empfangssignal 10 15 20 25 -0.5 0 0.5 1 (b) Variation von µLum ±1%. Zeit t/µs Skalierte Sensitivität ydssEL 10 15 20 25 -20 -10 0 10 20 (c) Skalierte Sensitivität bezüglich EL. Zeit t/µs Skalierte Sensitivität ydssµL 10 15 20 25 -4 -2 0 2 4 (d) Skalierte Sensitivität bezüglich µL. Abbildung 5.38: Skalierte Sensitivität ausgewählter elastischer Materialparameter auf das simulierte Empfangssignal. 5.6.2 Betrachtungen am viskoelastischen quasi-isotropen Material Ausgangspunkt für die Sensitivitätsanalyse des viskoelastischen Materials sei der im vorherigen Abschnitt betrachtete elastische Arbeitspunkt. Es wird der Fall der Quasi-Isotropie angenommen, entsprechend wird νL aus den anderen Materialparametern berechnet. Außerdem wird im Folgenden 5.6 Sensitivitätsanalyse 129 νLETELνTµL ycsspi/% 0 20 40 60 80 100 (a) Hexagonale Materialsymmetrie. ETELνTµLνL ycsspi/% 0 20 40 60 80 100 (b) Hexagonale Materialsymmetrie mit quasiisotroper Annahme. Abbildung 5.39: Informationsgehalt der simulierten diskreten Zeitsignal-Beobachtungen bezogen auf die elastischen Materialparameter ausgedrückt über die composite scaled sensitivity. µL nicht weiter betrachtet und jeweils nach Gl. (4.47) geschätzt. Neben dem Signalprototypen für das 750k-Schallwandler-Paar werden nun auch die vier weiteren Signalprototypen der übrigen Schallwandler-Paare genutzt, um der Frequenzabhängigkeit der Materialmatrix Rechnung zu tragen. Nach Gl. (4.54) werden für jede der vier Eigenbewegungen jeweils drei Parameter zur Beschreibung des fraktionalen Zener-Modells eingeführt. Da als Arbeitspunkt ein isotropes Material dient, ist entsprechend die Arbeitspunkt-Parametrisierung derart angepasst, dass die Scherkomponenten ( λ1 , λ2 , λ3 ) und die Dilatationskomponente ( λ4 ) unterschiedliche Verlustwinkel aufweisen, die Scherkomponenten untereinander jedoch denselben Verlustwinkel zugewiesen bekommen. Die exemplarische Parametrisierung ist in Tabelle 5.9 zusammengefasst. Auch in der folgenden Betrachtung der Sensitivität wird für die Bildung des zentralen Differenzenquotienten eine Parameteränderung von 1% genutzt. In Abb. 5.40(a) ist das Ergebnis der Sensitivitätsanalyse dargestellt, wobei aus Darstellungsgründen der Informationsgehalt bezüglich νi nicht aufgeführt wird. Es zeigt sich ein deutlich geringerer Informationsgehalt bezogen auf die Parameter des fraktionalen Zener-Modells (’viskoelastische’ Parameter) im Vergleich zum Informationsgehalt bezogen auf die ’elastischen’ Parameter. Wenig überraschend ist eine sehr geringe Sensitivität bezüglich der Viskoelastizität der Eigenbewegung 2 festzustellen, da diese dem Eigenwert λ−1 2=µL zugeordnet wird. Insgesamt zeigen alle drei Scher-Eigenbewegungen (1,2,3) den geringsten Informationsgehalt, sodass im Folgenden diese drei Eigenbewegungen bezüglich ihrer Verlustwinkel zusammengefasst werden. Eine mit dieser ModellVereinfachung wiederholte Betrachtung des Informationsgehalts ist in Abb. 5.40(b) dargestellt, wobei die Kennzeichnung G für Scher-Bewegungen und K für Dilatationsbewegungen steht. 130 5 Realisierung des inversen Messverfahrens Tabelle 5.9: Zusammenfassung der Parametrisierung der fraktionalen Zener-Modelle für die Scherund Dilatationskomponenten. Parameter τε/µsτσ/µsν Scherkomponente (Eigenwerte 1,2,3) 2,57 2,25 1 Dilatationskomponente (Eigenwert 4) 5,41 3,99 0,65 ETELνTτε4τσ4τε1τσ1τε3τσ3τε2τσ2 ycsspi/% 0 20 40 60 80 100 (a) Individuelle Betrachtung aller Eigenbewegungen. ETELνTτε,G τσ,G τε,K τσ,K νKνG ycsspi/% 0 20 40 60 80 100 (b) Reduzierte Betrachtung der Absorption auf Scherund Dilatationsbewegungen. Abbildung 5.40: Informationsgehalt der simulierten diskreten Zeitsignal-Beobachtungen bezogen auf die Materialparameter ausgedrückt über die composite scaled sensitivity. Reformulierung des viskoelastischen Materialmodells Auf Basis der in der obigen Analyse dargelegten Ergebnisse, wird das zu identifizierende Materialmodell vereinfacht zu ˜ C(iω) = "µTD1+µLD2+EL 1+νL D3,qi#1+ (iωτε,G)νG 1+ (iωτσ,G)νG +EL 1−2νL D4,qi 1+ (iωτε,K)νK 1+ (iωτσ,K)νK. (5.77) Anhand des vereinfachten Modells kann nun eine detailliertere Sensitivitätsanalyse erfolgen. Dabei zeigen sich folgende Auffälligkeiten bei der Betrachtung der skalierten Sensitivitäten, siehe Abb. 5.41: •τε,G , τσ,G , νG : Höhere Sensitivität auf den T-Puls 16 als auf den L-Puls 17 für Anregungen mit den Wandlern 750k und 1.0M. Aufgrund der starken Dämpfung des T-Pulses umgekehrtes Verhältnis der Sensitivitäten für die 1.5M-, 2.0Mund 2.5M-Wandler, siehe Abb. 5.41(a-b). 16 Signalanteile, welche nach dem Strahlenmodell von Reynolds mit den Anteilen der konvertierten Transversalwelle korrespondieren; Folge-Echos nach dem ersten Signal-Puls. 17 Signalanteile, welche nach dem Strahlenmodell von Reynolds mit den Anteilen der Longitudinalwelle korrespondieren; zeitlich zuerst eintreffender Signalanteil. 5.7 Inverses Problem 131 •τε,K , τσ,K , νK : Höhere Sensitivität auf den L-Puls als auf den T-Puls für alle Wandler-Paare, siehe Abb. 5.41(c-d). •νT , ET , EL : Höhere Sensitivität auf den L-Puls als auf den T-Puls für alle Wandler-Paare, siehe Abb. 5.41(e-h). • Sinkende Sensitivität aller untersuchten Parameter auf den Bereich des T-Pulses für steigende Frequenzen (Wandler höherer Mittenfrequenz), siehe Abb. 5.41(a-f). • Steigende Sensitivität aller untersuchten Parameter auf den Bereich des L-Pulses für steigende Frequenzen, siehe Abb. 5.41(a-f). 5.7 Inverses Problem In diesem Abschnitt wird das erarbeitete Verfahren zum Lösen des inversen Problems aufgezeigt. Zunächst werden Methoden zur Schätzung von Startwerten der Optimierungsvariablen (Messgrößen) beschrieben. Im Anschluss wird die Optimierungsstrategie diskutiert, welche die Auswahl des Optimierungsalgorithmus sowie die Formulierung der Kostenfunktion zusammenfasst. Zur Verdeutlichung der Zusammenhänge der im Folgenden diskutierten Signale sei definiert: xgen,j(t)7→ yp,j(t)7→ ymess,j(t) ysim,j(t) Generatorsignal Signalprototyp Empfangssignal Dabei sei der Index jnach Tabelle 5.10 definiert. Tabelle 5.10: Definition der Signal-Kennungen. jKennung Mittenfrequenz 1 750k 0,77 MHz 2 1.0M 1,15 MHz 3 1.5M 1,6 MHz 4 2.0M 2,1 MHz 5 2.5M 2,4 MHz 5.7.1 Startwertschätzverfahren Das allgemeine Vorgehen zur Startwertschätzung ist in zwei Schritte untergliedert. Zunächst werden für alle fünf Empfangssignale ymess,j(t) , welche sich durch unterschiedliche Mittenfrequenzen 132 5 Realisierung des inversen Messverfahrens Zeit t/µs Skalierte Sensitivität ydssτσ,G 10 15 20 25 -5 0 5 10 (a) Wandler-Paar 750k; Sensitivität auf τσ,G Zeit t/µs Skalierte Sensitivität ydssτσ,G 10 15 20 25 -15 -10 -5 0 5 10 15 (b) Wandler-Paar 1.5M; Sensitivität auf τσ,G Zeit t/µs Skalierte Sensitivität ydssτε,K 10 15 20 25 -6 -4 -2 0 2 4 6 (c) Wandler-Paar 750k; Sensitivität auf τε,K Zeit t/µs Skalierte Sensitivität ydssτε,K 10 15 20 25 -10 -5 0 5 10 (d) Wandler-Paar 1.5M; Sensitivität auf τε,K Zeit t/µs Skalierte Sensitivität ydssEL 10 15 20 25 -40 -20 0 20 40 (e) Wandler-Paar 750k; Sensitivität auf EL Zeit t/µs Skalierte Sensitivität ydssEL 10 15 20 25 -80 -60 -40 -20 0 20 40 60 (f) Wandler-Paar 1.5M; Sensitivität auf EL Zeit t/µs Skalierte Sensitivität ydssET 10 15 20 25 -60 -40 -20 0 20 40 60 (g) Wandler-Paar 750k; Sensitivität auf ET Zeit t/µs Skalierte Sensitivität ydssνT 10 15 20 25 -30 -20 -10 0 10 20 30 40 (h) Wandler-Paar 750k; Sensitivität auf νT Abbildung 5.41: Skalierte Sensitivität als Funktion der Zeit. 5.7 Inverses Problem 133 auszeichnen, Schätzwerte für die Longitudinalwellengeschwindigkeit, Transversalwellengeschwindigkeit, Longitudinalwellenabsorption und Transversalwellenabsorption berechnet und den jeweiligen Mittenfrequenzen zugeordnet. Nachfolgend werden diese frequenz-diskreten Kenngrößen verwendet, um ein isotropes viskoelastisches Materialmodell zu parametrisieren und somit einen Startparametervektor p(0)für das inverse Problem zu bilden. Signalanalyse im Verschiebezeitbereich Schätzung der Ausbreitungsgeschwindigkeit: In Abschnitt 5.2.3 wurde bereits das Strahlenmodell nach Reynolds zur Diskussion des Messeffekts und der Probekörperdimensionierung eingeführt. Hier soll nun das Strahlenmodell dazu dienen, geeignete Startwerte für die Longitudinalund Transversalwellengeschwindigkeit für ein Probenmaterial zu gewinnen. Als bekannt wird dazu die Geometrie der Probe vorausgesetzt. Nach Gl. (5.10) können Schätzwerte für cL und cT durch Auswertung der Eintreffzeitpunkte der Wellenpakete berechnet werden. In der Gleichung muss dazu der Außendurchmesser Da durch die Mantelstärke Dw ersetzt werden. Es gilt, den Eintreffzeitpunkt des ersten Wellenpakets (L-Puls) t0 und die Zeitdifferenz zwischen erstem und zweitem Wellenpaket (T-Puls) ∆t zu bestimmen, wobei zu beachten ist, dass die beiden Pulse zum einen dicht aufeinander folgen können und zum anderen der T-Puls sehr stark gedämpft sein kann. Methodisch wird das Problem durch Anwendung der analytischen Korrelation nach Gl. (B.21) gelöst, wobei zur Korrelation nicht das vom AWG ausgegebene Sendesignal xgen,j(t) , sondern der Signalprototyp yp,j(t) des jeweiligen Schallwandler-Paars genutzt wird. Der Signalprototyp beinhaltet bereits die gesamte Dynamik der Messkette, abgesehen von den Eigenschaften der Übertragungsstrecke (Probe), sodass realitätsnahe Schätzwerte mit bereits kompensierten systematischen Abweichungen berechnet werden können. Zur Berechnung des Sendeund Empfangsfrequenzgangs der zwei Schallwandler ist eine erste Schätzung bezüglich ihrer akustischen Last notwendig. Dazu wird die gravimetrisch ermittelte Dichte ρ=m/V der Probe und eine erste Schätzung der Longitudinalwellengeschwindigkeit ˘cL durch analytische Korrelation des Empfangssignals ymess,j(t) mit dem vom AWG ausgegebenen Sendesignal xgen,j(t) genutzt. Die akustische Last, welche den Schallwandlern zugeordnet wird, berechnet sich effektiv aus dem Flächenanteil, welcher durch die Probe mit der geschätzten akustischen Impedanz ˘ Zak belegt ist ( AProbe ) und dem Flächenanteil, welcher frei gegen Luft schwingt ( ALuft ). Die frei gegen Luft schwingende Fläche lässt sich dabei als Differenz zwischen der aktiven Schallwandlerfläche AWandler und der Probenfläche ausdrücken. Zeff,ak =ALuft430 Rayl +AProbe ˘ Zak mit ALuft =AWandler −AProbe (5.78) Nach Ausführen der analytischen Korrelation, wird der Betrag der Korrelationsfunktion nach lokalen Maxima durchsucht und zunächst die gefundenen Kandidaten entsprechend der Größe sortiert. Die größten zwei lokalen Maxima werden anschließend hinsichtlich ihrer Verschiebezeit sortiert.