The calculation of the dispersed melting behavior comes from the determination of the location of the first melting and the melting profile calculation together. The location of the first melting can be simplified to the location of the PFF. The carried out experiments validate the legitimacy of this simplification.

The physical, mathematical description of this melt process is assumed by the following requirements:

• The solid particles are assumed to be ideal not deformable spheres, which are evenly dispersed in the polymer melt.
• The particles are conceived as single particles. Changing effects between neighboring particles are ignored i.e. the surrounding area of a particle at-tunes to an undisturbed temperature field.
• The direction of the particle in the screw channel can be described by the central, cubic-faced sphere casing.
• The determination of the temperature rise in the solid area is assumed by the forced convection of the single particles. Temperature rise through heat radiation and –conduction will be neglected.
• The melting process begins at the location at which the solid particles are wetted by the melt. Melt formation in front of here is ignored.

For the training of a dispersed particle distribution at the location of the first filling, a minimum share of the already present melt $S_0$ is necessary. This is shown by:

$$S_0 = 1 - F_0 \tag{4.3}$$

Whereby $F_0$ represents the solid section. This solid section is identical to the bulk density and can be calculated using the assumptions of a central, cubic-faced particle formation.

During transportation within the solid conveying section, a rise in temperature of the solids particles takes place. At long solid conveying sections, the neglect of the extra heat content in this area could lead to mistakes in the calculation of the melting length. The average caloric particle temperature at the location of the first melting is consequently an input variable for the calculation of the melting profile.

If you consider the increase in temperature starting form a moving coordinate system (within a particle), the temperature rise poses an unsteady procedure. The essential heat transmitting mechanisms are heat conduction and convection.

Figure: Temperature trends within a spherically formed particle

In the area of unsteady heat transmittance, the heat increase starting from single particles [1, 2], or through solid beds [2, 3, 4], is considered. Since according to assumption 2, the melting of single particles should be assumed, it is also sensible to assume single particles in the solid conveying area. The figure shows diagrammatically the temperature trends inside a spherically formed particle. $\bar{T}$ is the average caloric temperature of the particle.

The coupling of the energy theorem and kinetics gives the differential calculus for the spherically symmetric temperature field.

$$\frac{\partial \theta}{\partial \tau} = \frac{1 \delta}{\xi^2 \delta \xi}\left(\xi^2 \frac{\delta \theta}{\delta \xi}\right) \tag{1}$$

whereby these standardization's:

$$\theta = \frac{T - T_0}{T_0 - T_U} ; \tau = \frac{at}{r_0^2} ; \xi = \frac{r}{r_0} \tag{2}$$

were introduced. For the case of a unique erratic temperature change from the starting temperature $T_0$ up to the environment temperature $T_U$, the solution of the differential calculus (2) is given by [2]. For sufficient lengths of time t in-side particles of finite expansion, similar temperature profiles can be expected. They are then described by the location function $f(x)$, which with extra time, experiences scaled reductions [2]:

$$\theta = g(\tau) \cdot f(\xi) \tag{3}$$

If one puts this product approach into the differential calculus so with the first order conditions one gets:

$$\theta(\xi, 0) = 1 \tag{4}$$

and the boundary conditions

$$\left(\frac{1}{Bi} \frac{\partial \theta}{\partial \xi} + \theta\right)_{\xi=1} = \begin{cases}1 \text{ für } \tau \leq 0\\0 \text{ für } \tau > 0\end{cases} \tag{5}$$

as a solution of the average caloric temperature of the particle:

$$\bar{\theta} = \sum_{i=1}^{\infty} c_i(m_i)D_i(m_i)e^{-m_i^2\tau} \tag{6}$$

with:

$$m_i = \left(1 - \frac{\alpha_L \cdot r_0}{\lambda}\right)\frac{\cos(m_i)}{\sin(m_i)} \tag{7}$$

$$c_i(m_i) = 2\frac{\sin(m_i) - m_i \cos(m_i)}{m_i - \sin(m_i)\cos(m_i)} \tag{8}$$

$$D_i(m_i) = 3\frac{\sin(m_i) - m_i \cos(m_i)}{m_i^3} \tag{9}$$

For fair calculation, it is sufficient to consider the first four elements of the sum function. The constants $m_i$ are the roots of the transcended eqn. (7). With the absolute Biot number

$$Bi = \frac{\alpha_L \cdot r_0}{\lambda} \tag{10}$$

The first four roots can be described through the following approximation equations:

$$m_1 = \frac{a}{\left(1 + \frac{b}{Bi}\right)^c} \tag{11}$$

$$m_{2-4} = a \cdot \tanh\left(\frac{\ln Bi}{\ln b} - c\right) + d \tag{12}$$

The constants a-d are stated in the table.

Table: Constants of the determination of the temperature function

The figure shows the value of the first four roots taken from [2] and the approximate values of those with the equations (11) and (12).

Figure: Roots for the calculation of the average caloric particle temperature

The dimensionless Bios number is dependent upon the external heat transfer coefficient $\alpha_L$, which can be determined with the help of the dimensionless Nusselt number.

$$Nu = \frac{\alpha_L \cdot \pi r_0}{\lambda_L} \tag{13}$$

The Nusselt number exposes itself from a laminate and turbulent section. For spherically formed single particles it is:

$$Nu = 2 + \sqrt{Nu_{lam}^2 + Nu_{tur}^2} \tag{14}$$

with:

$$Nu_{lam} = 0,664\sqrt{Re}\sqrt[3]{Pr} \tag{15}$$

$$Nu_{tur} = \frac{0,037Re^{0,8}Pr}{1 + 2,443Re^{-0,1}(Pr^{0,66} - 1)} \tag{16}$$

The dimensionless Reynolds number considers the air flow in the screw channel.

In the air direction, starting from a moving coordinate system in a particle, the as static assumed air moves relative to the particle within the screw channel.

It results in the air velocity being proportional to the conveying velocity of the particles.

$$Re = \frac{w\pi r_0}{v} \tag{17}$$

whereby $n_0$ is the screw speed and $t$ is the lead of the screw. The collection of solids in the screw channel reduces the flow effective channel cross section. This results in a rise in the flow velocity proportionate to the solid section $F_0$.

$$w = \frac{n_0 t}{F_0}$$

The Prandtl number contained in eqn. 15 and 16 is a ratio of material data. Using the here used consideration with air as a surrounding, round-flowing medium, the Prandtl number can be approximated by ignoring the pressure de-pendency using

$$\Pr(T)3,545 \cdot 10^{-7}T^2 - 1,309 \cdot 10^{-4}T + 0,7169 \tag{18}$$

Whereby the temperature has to be filled in with degrees Celsius.

Figure: Average caloric particle temperature as a function of location

The figure shows as an example for the material PE1810D and with an estimated particle diameter of 3mm the increase of the middle caloric particle temperature starting from the feeding temperature ($T_0$= 20 °C). For higher sur-rounding temperatures, higher average temperatures are determined. Assuming typical residence times in the solid conveying section of 1.0 – 1.5 seconds, it is clear to see that the increases in temperature cannot be ignored.

After the wetting of the solid particles by the melt, the particles can only be melted through heat conduction from the hot melt. The existing melt gets energy through dissipation and heat conduction.

The dispersed solid particles increase the yielded dissipation energy, since in this area of particles with ignored particle rotation, the shear gradient tends to zero. The consideration of these effects is conceived according to the figure, whereby the dispersion phase is seen as continual phase on the screw base.

Figure: Model of the consideration of the shear inflation

The flow effective channel depth is thus a function of the dimensionless solid content:

$$h(x,F) = (1 - F)h(x) \tag{1}$$

The calculation of the melt temperature development should be based on the following simplifications:

• Firstly, the screw channel is considered to be the shallow channel. The conditional change through the screw channel geometry in the shear heating of the melt compared to the shallow channel, are considered by an additionally inclusive corrective factor.
• The melt is wall-adhering
• The flow is stationary, incompressible and laminar slow moving ($c$ = $c_p$ = $c_v$)
• The flow characteristics of the melt follow the power law. The temperature dependency of the flow law's coefficients is described by Arrhenius' law: $\tau = K_{0T}e^{-\beta T}\dot{\gamma}^n$
• All material data, apart from the viscosity, are regarded as temperature in-variant.
• The dissipated energy per unit of volume is averaged over the channel cross-section and standardised on a reference temperature $T_j$.

With these assumptions, the general energy equation reduces to a description of the temperature field.

$$\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)} \tag{3}$$

The heat flow $q_y$ is represented by the Fourier heat conduction equation.

$$q_y = -\lambda \frac{\partial T}{\partial z} \tag{4}$$

When including the considered characteristic value, the described differential equation becomes:

$$\frac{\partial \theta}{\partial \zeta} = \frac{1}{Gz} \frac{\partial^2}{\partial \xi^2} + \frac{Br}{Gz} e^{-b T_z \theta} \tag{5}$$

Table: Dimensionless characteristics of the calculation of the melt temperature

For the reason of the exponential function, equation (5) is not solvable without further simplification. In order to make the problem solvable, the exponential function has been approximated by Potente [5] using a set of linear equations:

$$e^{-b(T-T_0)} \approx C_1 - T_z \beta(T - T_0) \tag{6}$$

The constants $C_1$ and $C_2$ are asserted for the specified areas by $\beta \cdot \Delta T$.

Figure: Approximation of the exponential function through a set of linear equations

In accordance with Potente, the stated constants are adjusted so bigger area values of $b^*DT$ can be allowed. The figure shows the approximation of the exponential function through equation 6. The table contains the constants $C_1$ and $C_2$ and their valid areas.

Table: Constants for the adaptation of the exponential function

For further analysis it should be assumed from this, that the temperature calculation is done in stages, and the temperature differences within the considered intervals $\Delta z$ are small. From this it follows that the expected temperature compensation processes are also small. This is sufficient if the change in heat flow in the channel direction by means of the start temperature profile is estimated, and is seen as an invariant opposite the z-coordinate. For every new interval, a new fixing takes place, because a new start temperature distribution occurs. The base temperature profile should suffice using the following conditions:

• No heat transfer to the screws occurs:

$$\frac{\partial \theta_0}{\partial \xi} = 0 \text{ für } \xi = 0 \tag{7}$$

• On the cylinder wall there is a given temperature $T_Z$:

$$\theta_0(\xi = 0) = \theta_z = \frac{T_z - T_0}{T_z} \tag{8}$$

• The average mass temperature of the melt at the beginning of the calculation section is known:

$$\int_0^1 \theta_0(\xi)d\xi = \theta_z = \frac{T_z - T_0}{T_z} \tag{9}$$

• An extreme should be permitted. This can be perceived by a large dissipation in the screw channel, or more specifically, in the melting zone.

A possible function that satisfies these conditions, is the superimposition of an exponential function with a second degree polynomial.

$$\theta_0(\xi) = 0.0154\theta_z(1 - e^{1-\xi})^8 - 0.1232\theta_z(\xi - 1)^2 - 0.04704\theta_z \tag{10}$$

The figure shows diagrammatically the starting temperature distribution for various cylinder wall temperatures according to equation (10).

Figure: Absolute base temperature distribution for various cylinder wall temperatures

With this base temperature profile as well as the acceptance that the change in heat flow in the deep channel direction can be estimated by means of the starting temperature profile:

$$\frac{\partial^2 \theta}{\partial \xi^2} \approx \frac{\partial^2 \theta_0}{\partial \xi^2} \tag{11}$$

you get the stated solution to the differential equation (5) in the table.

$$\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)$$

with $$\varepsilon = \frac{\theta_z}{Br C_2 \beta T_z}$$

Since only the channel height averaged mass temperature is of interest, this can be calculated by:

$$\bar{\theta} = \int_0^1 \theta_0(\xi)d\xi = \theta_z = \frac{T_z - T_0}{T_z} \tag{12}$$

The figure shows a comparison of the documented solutions in Table against the publicized approaches of Potente [4] and Ansahl [6]. For small Graetz numbers, the publicized approach of Ansahl tends to infinity, whilst both of the other curves converge on a different threshold value. At this point, the energy supplied from shearing equals the energy expelled from the heat conduction. For large Graetz numbers, subtract the deviation in the dimensionless temperature between the solutions.

Figure: Dimensionless mass temperature as a function of the Graetz number

The neglecting of the real channel geometry makes the analytical solution of the reduced energy equation ($\frac{\partial \theta}{\partial \zeta} = \frac{1}{Gz} \frac{\partial^2 \theta}{\partial \xi^2} + \frac{Br}{Gz} e^{-bT_z\theta}$) possible. For the reason of self-cleansing by the co-rotating twin screw extruder, the channel geometry varies greatly from the shallow channel. Through these simplifications, it can lead to big mistakes in the calculation of the average dissipated output per unit of volume $(\overline{\tau\dot{\gamma}})_j$, which moreover, still possess a marked dependency upon the local degree of filling. For the rectangular channel, the average dissipated output can be estimated by:

$$(\overline{\tau\dot{\gamma}})_j = \frac{K(T_j)v_0^{1+n}}{\bar{h}^{1+n}} \tag{1}$$

The changed averaged dissipated energy through the real channel geometry, as compared with the shallow channel, is taken into account by the correction factor $C_K$.

$$(\overline{\tau\dot{\gamma}})_j = C_K \frac{K(T_j)v_0^{1+n}}{\bar{h}^{1+n}} \tag{2}$$

For the determination of corrective factors, you have to approximate the real channel geometry using a step function.

Every interval i of the step function should possess, like the flat channel, a linear velocity distribution. If one lets the number of intervals tend towards infinity and the interval widths tend towards zero, it results in the average dissipated output:

$$(\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 \tag{3}$$

$x_f$ represents the position of the flow-front in the partially filled channel sections, assuming an ideal perpendicular flow-front. With the average effective channel depth

$$\bar{h}(x_f) = \frac{1}{x_f + \frac{b_{max}}{2}} \int_{-\frac{b_{max}}{2}}^{x_f} h(x)dx \tag{4}$$

solving for $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 \tag{5}$$

To avoid a null division for location $x=-b_{max}/2$ in equation (5), for the analogue to the output calculation, the practical channel profile should be used. Because of the segment definition of the channel profile, the determination of the intervals proves itself difficult, which is the reason for the integration being processed numerically. Equation (4) describes the filling degree dependent aver-age effective channel depth that has to be used for the plane channel model.

For the physical mathematical description of the melting of single particles in a polymer melt, according to prerequisite 2, changing effects between adjacent particles should be neglected. The energy equation in sphere coordination is reduced with the consideration of stationary relationships and natural heat at constant solid data to:

$$\rho c \frac{\partial T}{\partial t} = -\frac{1}{r^2} \frac{\partial}{\partial r}(r^2 \dot{q}_r) \tag{1}$$

To obtain the description of heat conductivity, Fourier's differential equation has to be used:

$$\dot{q}_r = -\lambda \frac{\partial T}{\partial r} \tag{2}$$

Put equation 2 into 1, under the assumption of constant material value to get:

$$\rho c \frac{\partial T}{\partial t} = \frac{\lambda}{r^2}\left(2r \frac{\partial T}{\partial t} + r^2 \frac{\partial^2 T}{\partial r^2}\right) \tag{3}$$

Figure: Sphere co-ordination on a solid particle

From the mass balance on the sphere face (see figure) it results in a change in the mass of the sphere per unit of time (Equation 4) equal to the change of the mass of the melt per unit of time (Equation (5).

Solid: $$\frac{\partial m_f}{\partial t} = -4\pi r_G^2 \rho_f \frac{\partial r_G}{\partial t} \tag{4}$$

Melt: $$\frac{\partial m_s}{\partial t} = -4\pi r^2 \rho_s \frac{\partial r}{\partial t} \tag{5}$$

By equating (4) and (5) one gets:

$$\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} \tag{6}$$

Solve equation (6) with respect to $\partial t$ and include the average flow velocity in the channel using $\left(\bar{v} = \frac{\partial z}{\partial t}\right)$ to get:

$$\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 \tag{7}$$

Now put equation (7) through (3) with:

$$a_s = \frac{\lambda_s}{\rho_s c_p} \tag{8}$$

With the definition of the constants:

$$\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 \tag{9}$$

From equation (8):

$$A = \frac{1}{a_s} r_G^2 \frac{\rho_f}{\rho_s} \bar{v} \frac{\partial r_G}{\partial z} \tag{10}$$

The double integration of

$$\frac{\partial^2 T}{\partial r^2} + \left[\frac{A}{r^2} + \frac{2}{r}\right] \frac{\partial T}{\partial r} = 0 \tag{11}$$

results in the following equation:

$$T(r) = \frac{C_1}{A}e^{-\frac{A}{r}} + C_2 \tag{12}$$

Whereby C1 and C2 are the integration constants. With the following boundary conditions:

$$T(r = \infty) = T_m \tag{13}$$

$$T(r = r_G) = T_{fl} \tag{14}$$

It results in the integration constants:

$$C_1 = A \frac{T_m - T_{fl}}{\exp\left(\frac{A}{r_G}\right) - 1} \tag{15}$$

$$C_2 = T_m + \frac{T_m - T_{fl}}{\exp\left(\frac{A}{r_G}\right) - 1} \tag{16}$$

The solution of the differential equation results in:

$$\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)} \tag{17}$$

Which is equilibrium of heat flows.

On the interface of the sphere, every time period $t \geq t_0$ has to be applied:

$\dot{q}_f|r_G = \dot{q}_s|r_G \tag{17}$

Wherein $\dot{q}_s$ is the melt-side heat-flow on the interface, and $\dot{q}_f$ is the heat-flow in the sphere on the interface.

For the heat-flow on the interface in the melt, with equation (16) the following is valid:

$\dot{q}_s = -\lambda \frac{\partial T}{\partial r} = \frac{\lambda(T_m - T_{fl})}{1 - \exp\left(\frac{A}{r_G}\right)}\left(-\frac{A}{r^2}\exp\left(\frac{A}{r}\right)\right) \tag{18}$

At the position r = $r_G$, for heat-flow the following is valid:

$\dot{q}_{s|r_G} = \lambda(T_m - T_{fl})\frac{A}{r_G^2} \frac{\exp\left(\frac{A}{r_G}\right)}{1 - \exp\left(\frac{A}{r_G}\right)} \tag{19}$

Then extend equation (20) with $\frac{\exp(-A/r_G)}{\exp(-A/r_G)}$ solving for heat flow results in:

$\dot{q}_{s|r_G} = \lambda(T_m - T_{fl})\frac{A}{r_G^2} \frac{1}{\exp\left(\frac{-A}{r_G}\right) - 1} \tag{20}$

Now substitute the constant A.

$A = \frac{1}{a_s} r_G^2 \frac{\rho_f}{\rho_s} \bar{v} \frac{\partial r_G}{\partial z} \tag{21}$

With the constant A':

$A' = \frac{1}{a_s} \frac{\rho_f}{\rho_s} \bar{v} \tag{22}$

So the heat-flow on the melt-side of the interface results in:

$\dot{q}_{s|r_G} = \lambda(T_m - T_{fl})A' \frac{\partial r_G}{\partial z} \frac{1}{\exp\left(-A'r_G \frac{\partial r_G}{\partial z}\right) - 1} \tag{23}$

For the heat-flow in the solid on the interface is valid:

$\dot{q}_{f|r=r_G} = \rho_f \bar{v}\Delta h \frac{\partial r_G}{\partial z} \tag{24}$

Then put equations (23) and (24) into equation 17 to get:

$\rho_f \bar{v}\Delta h = \frac{\lambda(T_m - T_{fl})A'}{\exp\left(-A'r_G \frac{\partial r_G}{\partial z}\right) - 1} \tag{25}$

Solve this equation with respect to $\frac{\partial r_G}{\partial z}$ to get:

$r_G \frac{\partial r_G}{\partial z} = -\frac{1}{A'}\ln\left[1 + \frac{\lambda A'(T_m - T_{fl})}{\rho_f \bar{v}\Delta h}\right] \tag{26}$

For the changing of the radius $r_0$ on a length $\Delta z$, through the integration of equation (26), the following is valid:

$r_{G,i+1} = \sqrt{r_{G,i}^2 - \frac{2}{A'}\ln\left[1 + \frac{\lambda A''(T_m - T_{fl})}{\Delta h}\right]\Delta z} \tag{27}$

include the constant A' in (27) to get:

$r_{G,i+1} = \sqrt{r_{G,i}^2 - \frac{2\lambda}{\rho_f \bar{v}c_s}\ln\left[1 + \frac{c_s(T_m - T_{fl})}{\Delta h}\right]\Delta z} \tag{28}$

$r_{0,i}$ is the radius at the beginning of the considered channel section. Equation (28) makes the calculation of the particle radius along the screw direction possible.

The solid section gives:

$F(r_G) = \frac{N_{p,ges} \frac{4}{3}\pi r_0^3}{\Delta z A_{channel}} \tag{29}$

Figure: Comparison of calculated and measured melt trends

The figure shows diagrammatically the comparison between experimentally determined (Symbols) and theoretically determined (linear) solid sections. The principal relation is accurately portrayed through the model. The solid section at the location of the first melt tallies with the experiment. What is notable is the calculated melt length is also accurately portrayed.