Dies ist eine alte Version des Dokuments!
Aufschmelzmodell für dispers verteilte Füllstoffe
Schmelzmodell für disperse Feststoffe
Die Berechnung des dispersen Schmelzverhaltens ergibt sich aus der Bestimmung des Ortes des ersten Schmelzens und der Berechnung des Schmelzprofils zusammen. Der Ort des ersten Schmelzens kann vereinfacht zum Ort des PFF werden. Die durchgeführten Experimente validieren die Berechtigung dieser Vereinfachung. Die physikalische, mathematische Beschreibung dieses Schmelzprozesses wird durch die folgenden Anforderungen angenommen:
- Die festen Partikel werden als ideale, nicht verformbare Kugeln angenommen, die gleichmäßig in der Polymerschmelze dispergiert sind.
- Die Partikel werden als Einzelpartikel betrachtet. Wechselwirkungen zwischen benachbarten Partikeln werden ignoriert, d. h. die Umgebung eines Partikels stimmt sich auf ein ungestörtes Temperaturfeld ab.
- Die Richtung des Partikels im Schneckenkanal kann durch die zentrale, kubisch-flächige Kugelhülle beschrieben werden.
- Die Bestimmung des Temperaturanstiegs im festen Bereich wird durch die erzwungene Konvektion der Einzelpartikel angenommen. Temperaturanstieg durch Wärmestrahlung und -leitung wird vernachlässigt.
- Der Schmelzprozess beginnt an dem Ort, an dem die festen Partikel von der Schmelze benetzt werden. Schmelzbildung vor diesem Ort wird ignoriert.
Für die Bildung einer dispersen Partikelverteilung am Ort der ersten Füllung ist ein minimaler Anteil der bereits vorhandenen Schmelze S0 notwendig. Dies wird dargestellt durch: $$S_0 = 1−F_0 \tag{4.3}$$ Wobei $F_0$ den Feststoffanteil darstellt. Dieser Feststoffanteil ist identisch mit der Schüttdichte und kann unter der Annahme einer zentralen, kubisch-flächigen Partikelbildung berechnet werden.
Feststofftemperatur am Ort des Aufschmelzbeginns
Während des Transports innerhalb der Feststoffförderzone findet eine Temperaturerhöhung der Feststoffpartikel statt. Bei langen Feststoffförderzonen kann die Vernachlässigung der Enthalpiezunahme in diesem Bereich zu Fehlern bei der Berechnung der Aufschmelzlänge führen. Die mittlere kalorische Partikeltemperatur am Ort des Aufschmelzbeginns ist somit eine Eingangsgröße für die Berechnung des Aufschmelzverlaufs.
Betrachtet man die Temperaturerhöhung ausgehend von einem bewegten Koordinatensystem (innerhalb eines Partikels), so stellt die Temperaturerhöhung einen instationären Vorgang dar. Die wesentlichen Wärmeübertragungsmechanismen sind Wärmeleitung und Konvektion.
Abbildung: Schematische Darstellung des Temperaturverlaufs im Innern eines kugelförmigen Partikels
Auf dem Gebiet der instationären Wärmeübertragung wird die Temperaturerhöhung ausgehend von Einzelpartikeln [1], [2] oder durchströmten Festbetten [2] [3] [4] betrachtet. Da nach Voraussetzung 2 für das Aufschmelzen von Einzelpartikeln ausgegangen werden soll, ist es sinnvoll, auch im Feststoffförderbereich von Einzelpartikeln auszugehen. Das Bild zeigt schematisch den Temperaturverlauf im Innern eines Partikels in Kugelform. $\bar{T}$ ist die gesuchte mittlere kalorische Temperatur der Partikel.
Die Kopplung von Energiesatz und Kinetik liefert die zu lösende Differentialgleichung für das kugelsymmetrische Temperaturfeld,
$$\frac{\partial\theta}{\partial\tau} = \frac{1}{\xi^2}\frac{\partial}{\partial\xi}\left(\xi^2\frac{\partial\theta}{\partial\xi}\right)$$
wobei die Normierungen
$$\theta = \frac{T - T_U}{T_0 - T_U}; \tau = \frac{at}{r_0^2}; \xi = \frac{r}{r_0}$$
eingeführt wurden. Für den Fall einer einmaligen sprunghaften Temperaturänderung von der Anfangstemperatur $T_0$ auf die Umgebungstemperatur $T_U$ ist in [2] die Lösung der Dgl. angegeben. Für hinreichend lange Zeiten $\tau$ im Innern der Partikel endlicher Ausdehnung sind die Temperaturprofile einander ähnlich. Sie lassen sich dann durch eine Ortsfunktion $f(\xi)$ beschreiben, die mit zunehmender Zeit maßstäblich verkleinert wird [2]:
$$\theta = g(\tau) \cdot f(\xi)$$
Führt man diesen Produktansatz in die Differentialgleichung ein, so erhält man mit der Anfangsbedingung
$$\theta(\xi,0) = 1$$
und der Randbedingung
$$\left(\frac{1\ \partial\theta}{\partial i\ \partial\xi} + \theta\right)_{\xi=1} = \begin{cases} 1 & \text{für } \tau \leq 0 \\ 0 & \text{für } \tau > 0 \end{cases}$$
als Lösung die kalorisch mittlere Temperatur der Partikel
$$\bar{\theta} = \sum_{i=1}^{\infty} c_i(m_i)D_i(m_i)e^{-m_i^2\tau}$$
mit
$$m_i = \left(1 - \frac{\alpha_L \cdot r_0}{\lambda}\right)\frac{\cos(m_i)}{\sin(m_i)}$$
$$c_i(m_i) = 2\frac{\sin(m_i) - m_i \cos(m_i)}{m_i - \sin(m_i) \cos(m_i)}$$
$$D_i(m_i) = 3\frac{\sin(m_i) - m_i \cos(m_i)}{m_i^3}$$
Für praxisgerechte Berechnungen ist es ausreichend, die ersten vier Glieder der Summenfunktion zu berücksichtigen. Die Konstanten $m_i$ sind die Nullstellen der in der ersten Gleichung angegebenen transzendenten Funktion. Mit der dimensionslosen Biot-Kennzahl
$$Bi = \frac{\alpha_L \cdot r_0}{\lambda}$$
können die ersten vier Nullstellen durch folgende Approximationsgleichungen beschrieben werden:
$$m_i = \frac{a}{\left(1 + \frac{b}{Bi}\right)^c}$$
$$m_{2-4} = a \cdot \tanh\left(\frac{\ln Bi}{\ln b} - c\right) + d$$
Die Konstanten a-d sind in der Tabelle angegeben.
| $a$ | $b$ | $c$ | $d$ | |
|---|---|---|---|---|
| $m_1$ | 3,140 | 10/3 | 1/2 | - |
| $m_2$ | 0,859 | 5 | 0,85 | 5,35 |
| $m_3$ | 0,875 | 7 | 1,1 | 8,6 |
| $m_4$ | 0,831 | 7 | 1,2 | 11,735 |
Tabelle: Konstanten zur Bestimmung der Temperaturfunktion
Folgendes Bild zeigt die Werte der ersten 4 Nullstellen aus [2] und die mit obenstehenden Gleichungen approximierten Werte.
Abbildung: Nullstellen zur Berechnung der mittleren kalorischen Partikeltemperatur
Die dimensionslose Biot-Zahl ist abhängig vom äußeren Wärmeübergangskoeffizienten $\alpha_L$, der sich mit Hilfe der dimensionslosen Nusselt-Zahl berechnen lässt.
$$Nu = \frac{\alpha_L \cdot \pi r_0}{\lambda_L}$$
Die Nusselt-Zahl setzt sich zusammen aus einem laminaren und einem turbulenten Anteil. Für kugelförmige Einzelpartikel gilt:
$$Nu = 2 + \sqrt{Nu_{lam}^2 + Nu_{tur}^2}$$
mit
$$Nu_{lam} = 0,664\sqrt{Re}\sqrt[3]{Pr}$$
$$Nu_{tur} = \frac{0,037Re^{0.8}Pr}{1 + 2,443Re^{-0.1}(Pr^{0.66} - 1)}$$
Die dimensionslose Reynolds-Zahl berücksichtigt die Luftströmung im Schneckenkanal.
$$Re = \frac{w\pi r_0}{\nu}$$
Betrachtet man die Luftgeschwindigkeit ausgehend von einem bewegten Koordinatensystem in einem Partikel, so bewegt sich die als ruhend im Schneckenkanal angenommene Luft relativ zum Partikel.
Es ergibt sich die Luftgeschwindigkeit proportional der Fördergeschwindigkeit des Granulates,
$$w = \frac{n_0 t}{F_0}$$
wobei $n_0$ die Schneckenantriebsdrehzahl und t die Gangsteigung der Schnecken ist. Die Feststoffansammlung im Schneckenkanal reduziert den strömungswirksamen Kanalquerschnitt. Hieraus resultiert eine Erhöhung der Strömungsgeschwindigkeit entsprechend des Feststoffanteils $F_0$.
Die in den Gleichungen (4.18) und (4.19) enthaltene Prandtl-Zahl Pr ist ein Stoffwert-Verhältnis. Bei der hier durchgeführten Betrachtung mit Luft als umströmendes Medium kann die Prandtl-Zahl bei Vernachlässigung der Druckabhängigkeit mit
$$Pr(T)3,545 \cdot 10^{-7}T^2 - 1,309 \cdot 10^{-4}T + 0,7169$$
approximiert werden, wobei die Temperatur in Grad Celsius einzusetzen ist.
Abbildung: Mittlere kalorische Partikeltemperatur als Funktion der Verweilzeit
Das Bild zeigt beispielhaft für das Material PE 1810 D und einem angenommenen Partikeldurchmesser von 3 mm die Erhöhung der mittleren kalorischen Partikeltemperatur ausgehend von der Zugabetemperatur ($T_0$=20°C). Für höhere Umgebungstemperaturen werden höhere mittlere Temperaturen ermittelt. Legt man typische Verweilzeiten in der Feststoffförderzone von 1.0 bis 1.5 Sekunden zugrunde, so erkennt man, dass die Temperaturerhöhung nicht vernachlässigt werden darf.
Schmelzetemperaturentwicklung
Nach der Benetzung der Feststoffpartikel mit Schmelze können diese nur durch Wärmeleitung von der heißen Schmelze aufgeschmolzen werden. Der vorhandenen Schmelze wird Energie durch Dissipation und Wärmeleitung zugeführt.
Die dispergierten Feststoffpartikel erhöhen die eingebrachte Dissipationsenergie, da im Bereich der Partikel bei vernachlässigter Partikelrotation der Schergradient zu Null wird. Zur Berücksichtigung dieses Effekts wird gemäß der Abbildung die disperse Phase als kontinuierliche Phase am Schneckengrund aufgefasst.
Abbildung: Modell zur Berücksichtigung der Scherüberhöhung
Die strömungswirksame Kanaltiefe ist somit eine Funktion des dimensionslosen Feststoffgehalts:
$$h(x,F) = (1 - F)h(x)$$
Für die Berechnung der Schmelzetemperaturentwicklung soll von folgenden Vereinfachungen ausgegangen werden:
- Der Schneckenkanal wird zunächst als Flachkanal betrachtet. Die durch die Schneckenkanalgeometrie bedingte Änderung in der Schererwärmung der Schmelze gegenüber dem Flachkanal wird durch einen nachträglich eingeführten Korrekturfaktor berücksichtigt.
- Die Schmelze ist wandhaftend.
- Die Strömung ist stationär, inkompressibel und laminar schleichend ($c$ = $c_p$ = $c_v$).
- Das Fließverhalten der Schmelze folgt dem Potenzgesetz. Die Temperaturabhängigkeit des Fließgesetzkoeffizienten wird mit dem Arrheniusansatz beschrieben: $\tau = K_{0T}e^{-\beta T}\dot{\gamma}^n$
- Alle Stoffwerte außer der Viskosität werden als temperaturinvariant betrachtet.
- Die dissipierte Leistung pro Volumeneinheit wird über den Kanalquerschnitt gemittelt und auf eine Bezugstemperatur $T_j$ normiert.
Mit diesen Annahmen reduziert sich die allgemeine Energiegleichung zur Beschreibung des Temperaturfeldes:
$$\frac{\partial T}{\partial z} = -\frac{1}{\rho_m c_p \bar{v}_z} \frac{\partial q_y}{\partial y} + \frac{(\tau\dot{\gamma})_j}{\rho_m c_p \bar{v}_z} e^{-\beta(T-T_j)}$$
Für den Wärmestrom gilt das Fourier'sche Wärmeleitgesetz:
$$q_y = -\lambda \frac{\partial T}{\partial z}$$
Die beschreibenden Dgl. lässt sich mit den in der Tabelle angegebenen Kennzahlen in eine dimensionslosen Darstellung überführen.
| $\frac{\partial \theta}{\partial \zeta} = \frac{1}{Gz} \frac{\partial^2 }{\partial \xi^2} + \frac{Br}{Gz} e^{-\beta T_Z \theta}$ | |
| $\theta = \frac{T - T_j}{T_z}$ | $\theta_0 = \frac{T_0 - T_j}{T_z}$ |
| $\xi= \frac{y}{h}$ | $\xi= \frac{z}{\Delta z}$ |
| $Br = \frac{\overline{\tau \cdot \dot{y}} \, \bar{h}^2}{\lambda T_z}$ | $Gz = \frac{c_p \rho \overline{h} \, }{\lambda \overline{b} \Delta z} \dot{V}_z$ |
Tabelle: Dimensionslose Kenngrößen für die Berechnung der Schmelzetemperatur
Die Gleichung ist aufgrund der Exponentialfunktion nicht ohne weitere Vereinfachungen lösbar. Um das Problem lösbar zu machen, ist die Exponentialfunktion von Potente [5] durch ein Polygonzug angenähert worden:
$$e^{-\beta(T-T_0)} \approx C_1 - T_2\beta(T - T_0)$$
Die Konstanten $C_1$ und $C_2$ sind für bestimmte Bereiche von $\beta \cdot \Delta T$ fest vorgegeben.
Abbildung: Annäherung der Exponentialfunktion durch einen Polygonzug
In Anlehnung an Potente werden die Konstanten hier neu angepasst, um einen größeren Wertebereich von $\beta \cdot \Delta T$ zulassen zu können. Das Bild zeigt die Annäherung der Exponentialfunktion durch die Gleichung, nachstehende Tabelle enthält die Konstanten $C_1$ und $C_2$ und deren Gültigkeitsbereich.
| $\beta \cdot \Delta T$ | $c_1$ | $c_2$ |
|---|---|---|
| $-0.5 \cdot \Delta T \leq \beta \cdot \Delta T$ | 1.0289 | 1.0508 |
| $0.2 \leq \beta \cdot \Delta T \leq 0.8$ | 0.9419 | 0.6157 |
| $0.8 \leq \beta \cdot \Delta T < 1.5$ | 0.7322 | 0.3536 |
Tabelle: Konstanten zur Anpassung der Exponentialfunktion
Für die weitere Analyse soll davon ausgegangen werden, dass die Temperaturberechnung abschnittsweise erfolgt und die Temperaturdifferenzen innerhalb des betrachteten Intervalls $\Delta z$ klein sind. Hieraus lässt sich folgern, dass die zu erwartenden Temperaturausgleichsvorgänge ebenfalls klein sind. Es ist daher ausreichend, wenn die Änderung des Wärmestroms in Kanaltiefenrichtung anhand des Anfangstemperaturprofils abgeschätzt wird und als invariant gegenüber der z-Koordinate angesehen wird. Für jedes neue Intervall findet eine neue Festlegung statt, weil sich eine neue Anfangstemperaturverteilung ergibt.
Das Anfangstemperaturprofil soll dabei folgenden Bedingungen genügen:
Es erfolgt kein Wärmeübergang zur Schnecke:
$$\frac{\partial\theta_0}{\partial\xi} = 0$$ für $\xi = 0$
An der Zylinderwand liegt eine vorgegebene Temperatur $T_Z$ vor:
$$\theta_0(\xi = 0) = \theta_z = \frac{T_z - T_0}{T_z}$$
Die mittlere Massetemperatur der Schmelze zu Beginn des Berechnungsabschnittes ist bekannt:
$$\int_0^1 \theta_0(\xi)d\xi = \theta_s = \frac{T_s - T_0}{T_z}$$
Ein Extremum soll zugelassen sein. Dies ist bei großer Dissipation im Schneckenkanal z. B. in der Aufschmelzzone denkbar.
Eine mögliche Funktion, die diesen Bedingungen genügt, ist die Überlagerung einer Exponentialfunktion mit einen Polynom 2. Grades.
$$\theta_0(\xi) = 0.0154\theta_z(1 - e^{1-\xi})^8 - 0.1232\theta_z(\xi - 1)^2 - 0.04704\theta_z$$
Folgendes Bild zeigt beispielhaft die Anfangstemperaturverteilung gemäß der Gleichung für verschiedene dimensionslose Zylinderwandtemperaturen.
Abbildung: Dimensionslose Anfangstemperaturverteilung für verschiedene dimensionslose Zylinderwandtemperaturen
Mit diesem Anfangstemperaturprofil sowie der Annahme, dass die Änderung des Wärmestroms in Kanaltiefenrichtung anhand des Anfangstemperaturprofils abgeschätzt werden kann,
$$\frac{\partial^2\theta}{\partial\xi^2} \approx \frac{\partial^2\theta_0}{\partial\xi^2}$$
erhält man die in der Tabelle angegebene Lösung der Dgl.
$$\theta(\xi,\zeta) = \frac{C_1}{C_2\beta T_z} - 0,2464\varepsilon - 0,1232\varepsilon e^{1-\xi} + 1,7248\varepsilon e^{2(1-\xi)} - 7,7616\varepsilon e^{3(1-\xi)} + 17,2479\varepsilon e^{4(1-\xi)}$$
$$-21,5599\varepsilon e^{5(1-\xi)} + 15,5231\varepsilon e^{6(1-\xi)} - 6,0368\varepsilon e^{7(1-\xi)} + 0,9856\varepsilon e^{8(1-\xi)}$$
$$-e^{\frac{BrC_2\beta T_z\zeta}{Gz}}\left(\frac{C_1}{C_2\beta T_z} - 0,2464\varepsilon - 0,1232\varepsilon e^{1-\xi} + 1,7248\varepsilon e^{2(1-\xi)} - 7,7616\varepsilon e^{3(1-\xi)}\right.$$
$$+ 17,2479\varepsilon e^{4(1-\xi)} - 21,5599\varepsilon e^{5(1-\xi)} + 15,5231\varepsilon e^{6(1-\xi)} - 6,0368\varepsilon e^{7(1-\xi)}$$
$$+ 0,9856\varepsilon e^{8(1-\xi)} + 0,1548\theta_z + 0,1232\theta_z e^{1-\xi} - 0,4312\theta_z\left(e^{1-\xi}\right)^2$$
$$+ 0,8624\theta_z e^{3(1-\xi)} - 1,0780\theta_z e^{4(1-\xi)} + 0,8624\theta_z e^{5(1-\xi)} - 0,4312\theta_z e^{6(1-\xi)}$$
$$\left.+ 0,1232\theta_z e^{7(1-\xi)} - 0,0154\theta_z e^{8(1-\xi)} + 0,1232\theta_z\xi^2 - 0,2464\theta_z\xi\right)$$
mit
$$\varepsilon = \frac{\theta_s}{BrC_2\beta T_z}$$
Da in der Regel nur eine über die Kanalhöhe gemittelte Massetemperatur von Interesse ist, kann diese mit
$$\bar{\theta} = \frac{\bar{T} - T_0}{T_z} = \int_0^1 \theta(\xi,\zeta)d\xi$$
bestimmt werden.
Nachstehendes Bild zeigt einen Vergleich der in der Tabelle dokumentierten Lösung mit den von Potente [4] und Ansahl [6] publizierten Ansätzen. Für kleine Graetz-Zahlen strebt die von Ansahl publizierte Lösung ins Unendliche, während die beiden anderen Kurven auf unterschiedliche Grenzwerte zulaufen, bei denen die durch Scherung zugeführte Energie gleich der durch Wärmeleitung abgeführten Energie ist. Für große Graetz-Zahlen nehmen die Abweichungen in der dimensionslosen Temperatur zwischen den Lösungen ab.
Abbildung: Dimensionslose Massetemperatur als Funktion der Graetz-Zahl
Berücksichtigung der realen Kanalgeometrie
Die Vernachlässigung der realen Kanalgeometrie ermöglicht die analytische Lösung der reduzierten Energiegleichung. Aufgrund der Selbstreinigung von Gleichdrall-Doppelschneckenextrudern weicht jedoch die Kanalgeometrie stark von dem Flachkanal ab. Durch diese Vereinfachung kann es zu größeren Fehlern bei der Bestimmung der gemittelten dissipierten Leistung pro Volumeneinheit $\overline{(\tau\dot{\gamma})_j}$ kommen, die zudem noch eine ausgeprägte Abhängigkeit vom örtlichen Füllgrad besitzt.
Für den flächengleichen Rechteckkanal kann die gemittelte dissipierten Leistung mit $\overline{(\tau\dot{\gamma})_j} = \frac{K(T_j)v_0^{1+n}}{\bar{h}^{1+n}}$$
abgeschätzt werden. Die durch die reale Kanalgeometrie bedingte Änderung der gemittelten dissipierten Leistung gegenüber dem Flachkanal soll durch einen zusätzlich eingeführten Korrekturfaktor $C_K$ Berücksichtigung finden.
$$\overline{(\tau\dot{\gamma})_j} = C_K \frac{K(T_j)v_0^{1+n}}{\bar{h}^{1+n}}$$
Zur Ermittlung des Korrekturfaktors nähert man die reale Kanalgeometrie durch eine Treppenfunktion an. Jedes Intervall i der Treppenfunktion möge wie der Flachkanal eine lineare Geschwindigkeitsverteilung besitzen. Lässt man die Anzahl der Intervalle gegen unendlich und damit die Intervallbreite gegen Null streben, so ergibt sich für die gemittelte dissipierte Leistung:
$$(\overline{\tau\dot{\gamma}})_j = \frac{1}{x_f + \frac{b_{max}}{2}} \int_{-\frac{b_{max}}{2}}^{x_f} \frac{K(T_j)v_0^{1+n}}{h(x)^{1+n}} dx$$
$x_f$ ist dabei die Position der Fließfront in teilgefüllten Kanalabschnitten, wobei von einer idealisierten senkrechten Fließfront ausgegangen wird. Mit der mittleren wirksamen Kanaltiefe
$$\bar{h}(x_f) = \frac{1}{x_f + \frac{b_{max}}{2}} \int_{-\frac{b_{max}}{2}}^{x_f} h(x)dx$$
erhält man für $C_K$:
$$C_K = \frac{1}{x_f + \frac{b_{max}}{2}} \int_{-\frac{b_{max}}{2}}^{x_f} \left(\frac{\bar{h}(x_f)}{h(x)}\right)^{1+n} dx$$
Damit in der Gleichung eine Nulldivision für den Ort x=-$b_{max}$/2 vermieden wird, muss auch hier analog zur Leistungsberechnung das praktische Kanalprofil verwendet werden. Die Ermittlung des Integrals erweist sich aufgrund der abschnittsweisen Definition des Kanalprofils als schwierig, weshalb die Integration numerisch durchgeführt wurde. Die Gleichung beschreibt die füllgradabhängige mittlere wirksame Kanaltiefe, die für das ebene Rinnenmodell zu verwenden ist.
Berechnung der Feststoffbettabnahme entlang des Aufschmelzweges
Für die physikalisch mathematische Beschreibung des Aufschmelzens von Einzelpartikeln in einer Polymerschmelze sollen gemäß Voraussetzung 2 Wechselwirkungen zwischen benachbarten Partikeln vernachlässigt werden. Die Energiegleichung in Kugelkoordinaten reduziert sich bei Betrachtung stationärer Verhältnisse und reiner Wärmeleitung bei konstanten Stoffdaten:
$$\rho c \frac{\partial T}{\partial t} = -\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\dot{q}_r\right)$$
Zur Beschreibung der Wärmeleitung verwendet man die Fourier'sche Differentialgleichung
$$\dot{q}_r = -\lambda \frac{\partial T}{\partial r}$$
Setzt man die Gleichungen ineinander ein, so erhält man unter der Annahme konstanter Stoffwerte:
$$\rho c \frac{\partial T}{\partial t} = \frac{\lambda}{r^2}\left(2r\frac{\partial T}{\partial r} + r^2\frac{\partial^2 T}{\partial r^2}\right)$$
Abbildung: Kugelkoordinaten am Feststoffpartikel
Aus der Massebilanz an der Kugeloberfläche ergibt sich, dass die Änderung der Masse der Kugel pro Zeiteinheit gleich der Änderung der Masse der Schmelze pro Zeiteinheit ist.
Feststoff: $$\frac{\partial m_f}{\partial t} = -4\pi r_G^2 \rho_f \frac{\partial r_G}{\partial t}$$
Schmelze: $$\frac{\partial m_s}{\partial t} = -4\pi r^2 \rho_s \frac{\partial r}{\partial t}$$
Durch Gleichsetzen erhält man: $$\frac{\partial r}{\partial t} = -\left(\frac{r_G}{r}\right)^2 \frac{\rho_f}{\rho_s} \frac{\partial r_G}{\partial t} = -\left(\frac{r_G}{r}\right)^2 \frac{\rho_f}{\rho_s} \frac{\partial r_G}{\partial_z} \frac{\partial_z}{\partial t}$$
Löst man die Gleichung nach $\partial t$ auf und führt die mittlere Strömungsgeschwindigkeit im Kanal mit $\bar{v} = \frac{\partial z}{\partial t}$ ein, so erhält man:
$$\partial t = \left(-\left(\frac{r_G}{r}\right)^2 \frac{\rho_f}{\rho_s} \bar{v} \frac{\partial r_G}{\partial z}\right)^{-1} \partial r$$
durch Einsetzen erhält man:
$$\frac{\partial^2 T}{\partial r^2} + \left[\frac{1}{a_s}\left(\frac{r_G}{r}\right)^2 \frac{\rho_f}{\rho_s} \bar{v} \frac{\partial r_G}{\partial z} + \frac{2}{r}\right] \frac{\partial T}{\partial r} = 0$$
mit $$a_s = \frac{\lambda_s}{\rho_s c_p}$$
Mit der Definition der Konstanten:
$$A = \frac{1}{a_s}r_G^2\frac{\rho_f}{\rho_s}\bar{v}\frac{\partial r_G}{\partial z}$$
erhält man:
$$\frac{\partial^2 T}{\partial r^2} + \left[\frac{A}{r^2} + \frac{2}{r}\right]\frac{\partial T}{\partial r} = 0$$
Die zweimalige Integration führt zu Gleichung:
$$T(r) = \frac{C_1}{A}e^{-\frac{A}{r}} + C_2$$
wobei $C_1$ und $C_2$ die Integrationskonstanten sind. Mit den Randbedingungen:
$$T(r = \infty) = T_m$$ $$T(r = r_G) = T_{fl}$$
ergeben sich die Integrationskonstanten zu:
$$C_1 = A \frac{T_m - T_{fl}}{\exp\left(\frac{A}{r_G}\right) - 1}$$
$$C_2 = T_m + \frac{T_m - T_{fl}}{\exp\left(\frac{A}{r_G}\right) - 1}$$
Die Lösung der Differentialgleichung lautet:
$$\frac{T_m - T(r)}{T_m - T_{fl}} = \frac{1 - \exp\left(\frac{A}{r}\right)}{1 - \exp\left(\frac{A}{r_G}\right)}$$
Wärmestrombilanz
An der Grenzfläche der Kugel muss zu jeder Zeit $t \geq t_0$ gelten:
$$\dot{q}_f|_{r_G} = \dot{q}_s|_{r_G}$$
wobei $\dot{q}_s$ der Wärmestrom schmelzeseitig an der Grenzfläche und $\dot{q}_f$ der Wärmestrom in der Kugel an der Grenzfläche ist.
Für den Wärmestrom an der Grenzfläche in der Schmelze ergibt sich:
$$\dot{q}_s = -\lambda \frac{\partial T}{\partial r} = \frac{\lambda\left(T_m - T_{fl}\right)}{1 - \exp\left(\frac{A}{r_G}\right)}\left(-\frac{A}{r^2}\exp\left(\frac{A}{r}\right)\right)$$
An der Stelle r = $r_G$ ergibt sich damit für den Wärmestrom:
$$\dot{q}_s|_{r_G} = \lambda\left(T_m - T_{fl}\right)\frac{A}{r_G^2}\frac{\exp\left(\frac{A}{r_G}\right)}{1 - \exp\left(\frac{A}{r_G}\right)}$$
Erweitert man mit $\frac{\exp\left(-A/r_G\right)}{\exp\left(-A/r_G\right)}$, so erhält man für den Wärmestrom:
$$\dot{q}_s|_{r_G} = \lambda\left(T_m - T_{fl}\right)\frac{A}{r_G^2}\frac{1}{\exp\left(\frac{-A}{r_G}\right) - 1}$$
Ersetzt man nun die Konstante A
$$A = \frac{1}{a_s}r_G^2\frac{\rho_f}{\rho_s}\bar{v}\frac{\partial r_G}{\partial z}$$
durch die Konstante A':
$$A' = \frac{1}{a_s}\frac{\rho_f}{\rho_s}\bar{v}$$
so ergibt sich für den Wärmestrom an der Grenzfläche schmelzeseitig:
$$\dot{q}_{s|r_G} = \lambda\left(T_m - T_{fl}\right)A'\frac{\partial r_G}{\partial z}\frac{1}{\exp\left(-A'r_G\frac{\partial r_G}{\partial z}\right) - 1}$$
Für den Wärmestrom im Feststoff an der Grenzfläche gilt:
$$\dot{q}_{f|r=r_G} = \rho_f\bar{v}\Delta h\frac{\partial r_G}{\partial z}$$
Durch Einsetzen erhält man:
$$\rho_f\bar{v}\Delta h = \frac{\lambda\left(T_m - T_{fl}\right)A'}{\exp\left(-A'r_G\frac{\partial r_G}{\partial z}\right) - 1}$$
Löst man diese Gleichung nach $\frac{\partial r_G}{\partial z}$ auf, so erhält man:
$$r_G\frac{\partial r_G}{\partial z} = -\frac{1}{A'}\ln\left[1 + \frac{\lambda A'\left(T_m - T_{fl}\right)}{\rho_f\bar{v}\Delta h}\right]$$
Durch Integration erhält man für die Änderung des Radius' $r_0$ auf einer Länge $\Delta z$:
$$r_{G,i+1} = \sqrt{r_{G,i}^2 - \frac{2}{A'}\ln\left[1 + \frac{\lambda A''\left(T_m - T_{fl}\right)}{\Delta h}\right]\Delta z}$$
Setzt man in Gl. (4.66) A' ein, so ergibt sich:
$$r_{G,i+1} = \sqrt{r_{G,i}^2 - \frac{2\lambda}{\rho_f\bar{v}c_s}\ln\left[1 + \frac{c_s\left(T_m - T_{fl}\right)}{\Delta h}\right]\Delta z}$$
$r_{0,i}$ ist hierbei der Radius am Anfang des betrachteten Kanalabschnittes. Gleichung (4.67) ermöglicht somit die Berechnung des Partikelradius' entlang der Schneckenanordnung. Der Feststoffanteil ergibt sich zu:
$$F(r_G) = \frac{N_{p,ges} \frac{4}{3}\pi r_0^3}{\Delta z A_{channel}}$$
Abbildung: Vergleich berechneter und gemessener Aufschmelzverläufe
Das Bild zeigt beispielhaft den Vergleich zwischen experimentellen (Symbole) und theoretisch ermittelten (durchgezogene Linie) Feststoffanteilen. Das prinzipielle Verhalten wird durch das Modell richtig wiedergegeben. Der Feststoffanteil am Ort des Aufschmelzbeginns stimmt gut mit dem Experiment überein. Auch die berechneten Aufschmelzlängen werden gut wiedergegeben.