The feet material is held in underneath the opening by the screw-elements and is transported in the screws as a solid in the direction of the screw-tips. Through contact with the heated cylinder wall, the particles in the proximity of the cylinder wall begin to sinter and melt on until finally a broad standing connected molten film over the solid tailback in the intermeshing area is formed. This procedure is supported through in the practice conventionally carried-out reductions of the short pitch of the screw elements. These are found in front of the geometric plastification zone, through which an extra compression of the solids occurs and thus a better contact with the heated cylinder wall can be achieved. At the same time, granule particles lying on the underside of this layer are heated through convection.

Through the formation of the molten film, the friction conditions in the mesh area and in the solid tailback change behind the intermeshing area. Through this process, the forced conveying of solids in a intermeshing region breaks down, and material is forced through the nip region. As a result of the 3-dimensional speed profile in the intermeshing area, an intensive mixture of the materials occurs. If the proportion of ready-fused materials is sufficient, a dispersion from the existing granules is formed from the solid phase and from the existing viscous phase polymer melt.

If the melted share for a completely formed dispersion is not yet sufficient an agglomeration of granule particles occurs. These are then further fused through the dissipation of the already existing melt, and through heat conduction of the surrounding walls. Also the agglomerates are further more molten and trans-ported during their transportation through convection and through passing of further intermeshing areas in the already stated condition of dispersion.

After the achievement of this condition, the further melting of the not yet molten Polymer particles can only be successful through the heat conduction from the hot melt to the solid granule particles. This is because the granule particles no longer have direct contact with the heated cylinder walls, i.e. they are longer touching the heated screw surface.

A condition for the described melting behavior is a sufficiently long solids-conveying section. In practice, inserted screw configurations point out that with a short conveying section, combinations of plastic elements with subsequent dust particles, results in only a partial filling in the plastification zone.

If the available length of conveying section to the molten film over the solid tail-back does not suffice, so the forced conveying of the solids remains up to the location where the first filling with melt took place. Principally, if the solid that holds a significantly larger speed component in the axial direction than the melt, reaches the melt filled area, and melt can penetrate in the cavities of the granules, then the forced conveying will collapse. In front of the location of the first filling up with melt, through a sufficiently large mass throughput, originates one with a granule filled screw area, in which, the granules rotate around the screw like in the figure.

Figure: Melting behavior in a co-rotating twin screw extruder

The speed component in the axial direction is substantially lower in this area than in that of the solid-forced conveying area. Between that in the fully filled solid area emerging friction and the from the extrusion flow resulting force an equilibrium attunes.

The necessary pressure for wetting of particles can be worked out by:

$$p = \frac{2K(T_m)d_p}{a_f - d_p}\left(\frac{2F(1 + 2n)}{n(a_f - d_p)}\right)^n \tag{1}$$

When assuming quadratic gap areas.

Where $\bar{v}$ is the average flow velocity of the melt that arises through the evaluation of the continuity equation:

$$\bar{v} = \frac{\dot{m}}{A_{channel} \cdot \rho_s}\left(\frac{1}{1 - \frac{1}{F_0 \rho_m + (1 - F_0)\rho_m}}\right) \tag{2}$$

With equations (1) and (2) it is possible to calculate the length of the solid tailback.

In the solid filled area, through the heated cylinder walls, as well as with the hot screws, it can likewise lead to a molten film formation. A noticeable melt, however, sets in when solid particles enter the melt filled area. Here you will find that the particles are dispersed in the melt area. The further melting is then again successful through the heat conduction of the hot melt to the solid particles, whereby the particles diameter in the cross section are reduced.

After the achievement of this condition, the distribution of the not yet melted solid particles in the polymer melt is shown in the figure in the illustration of the thinly sliced samples. The cut direction runs parallel to the screw axis.

At the point of the first filling with melt, the solids rate is just over 50%. In this area, small air holes are still observable. What this points out is that not all cavities in this area have been completely filled with melt. In further progressions of melting these air holes should disappear. At the same time the size of the solid particles will reduce until they are non-existent. A preferred distribution of the solid particles in the polymer melt during the melting process is not perceivable.

For a global estimation of the extrusion process, the knowledge of the melting process is of a great importance. To calculate melting length and the solids bed profile, a modified TADMOR model is used. This estimates the location de-pendent molten film thickness on the cylinder wall amongst the inclusion of radial leakage flows. The solution of the melting process in the channel direction results from the energy balance in the molten film above the solid bed and the mass balance between the melt pool and the melt film, see figure.

Figure: Modified TADMOR model with location dependent melt film thickness

The equations for the calculation of the melt profile in the channel direction can be seen in the table.

Dimensionless solid base width from the energy balance in the melt film: $$y = \frac{X}{b} = \frac{\frac{k_1 \rho_s v_{0z} \Delta h}{\lambda_z(T_z - T_{Fl})b}\left(\frac{\delta}{1 + c}\right)^2}{2 + \frac{k_2 K(T_{Fl})v_{rel}^{1+n}}{\lambda_z(T_z - T_{Fl})}\left(\frac{\delta}{1 + c}\right)^{1-n}}\left(1 - \frac{s_R}{\delta}\right)\tag{3}$$

Dimensionless melt film width: $$\psi = \frac{\delta}{\delta_0}\psi^* = \frac{\psi - \psi_S}{1 - \psi_S} = y^c\tag{4}$$

Constant gradient for the melt profile: $$c = \frac{\lg\left(\frac{\psi_1 - \psi_S}{\psi_2 - \psi_S}\right)}{\lg\left(\frac{y_1}{y_2}\right)}\tag{5}$$

Function of the location dependent melt film width: $$\delta_0 = (\delta_1 - s_R)y_1^{-c} + s_R\tag{6}$$

Constants for the melt calculation: $$k_1 = 2\left(\frac{1}{1 - e^A} + \frac{1}{A}\right) k_2 = \frac{2}{A^2}\left(\frac{A}{e^A - 1}\right)^{1+n}(e^A - A - 1)\tag{7}$$ $$A = \frac{\beta}{n}(T_z - T_{Fl})\tag{8}$$

Function of the calculation of the dimensionless solid base width in the channel direction: $$y = \frac{X}{b} = [1 - (1 - c)(1 - \psi_S)\pi_1 \zeta]^{\frac{1}{1-c}}\tag{9}$$ $$\pi_1 = \frac{\rho_s k_1 \delta_0 v_0 D_S}{2\dot{m}}\zeta = \frac{z}{D_S} = \frac{L}{D_S \sin(\varphi_S)}\tag{10}$$

The location where the melting process starts can be simplified to the point of first filling (PFF). These assumptions were validated by various arrays of experiments.

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{11}$$

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 [Mar84], [NN88] or through solid beds [NN88], [TM87], [Sch84], 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{12}$$

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{13}$$

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 13 is given by [NN88]. 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 [NN88]:

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

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

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

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{16}$$

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

with:

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

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

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

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. 18. With the absolute Biot number

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

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

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

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

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 [NN88] and the approximate values of those with the equations 22 and 23.

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{24}$$

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{25}$$

with:

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

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

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{28}$$

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

The Prandtl number contained in eqn. 26 and 27 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{30}$$

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{31}$$

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{32}$$

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

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

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

Table: Dimensionless characteristics of the calculation of the melt temperature

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

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

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{35}$$

• 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{36}$$

• 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{37}$$

• 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{38}$$

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

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{39}$$

you get the stated solution to the differential equation 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)\tag{40}$$

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

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{42}$$

The figure shows a comparison of the documented solutions in Table against the publicized approaches of Potente [Pot90] and Ansahl [Ans93]. 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{43}$$

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{44}$$

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{45}$$

$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{46}$$

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{47}$$

To avoid a division by zero in the equation at the point $x=-b_{max}/2$, the practical channel profile must be used here as well, in the same way as for the power calculation. Determining the integral proves difficult due to the segmented definition of the channel profile, which is why the integration was performed numerically. The equation describes the flow depth-dependent mean effective channel depth, which is to be used for the flat 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{48}$$

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{49}$$

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{50}$$

Figure: Sphere co-ordination on a solid particle

The mass balance at the surface of the sphere shows that the change in the mass of the sphere per unit of time is equal to the change in the mass of the melt per unit of time.

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

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

Setting the two equal gives:

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

Solve equation 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{54}$$

Substituting gives:

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

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{56}$$

From equation 55:

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

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{58}$$

results in the following equation:

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

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

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

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

It results in the integration constants:

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

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

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{64}$$

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{65}$

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.

The heat flux at the interface within the melt is given by:

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

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{67}$

Then extend equation (68) 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{68}$

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{69}$

With the constant A':

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

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{71}$

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{72}$

Substituting gives:

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

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{74}$

For the changing of the radius $r_0$ on a length $\Delta z$, through the integration of equation (74), 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{75}$

include the constant A' in (75) 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{76}$

$r_{0,i}$ is the radius at the beginning of the considered channel section. Equation 76 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{77}$

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.

The modified disperse melting model introduced in [Thü08] is basically build on the theoretical remarks on disperse melting according to Melisch [Mel98], [PM96]. At this point references are made to the preliminary considerations introduced in the melting model for disperse fillers.

In the following the assumptions and the boundary conditions introduced are used to find a detailed solution.

The influence of the finite channel dimension is especially pronounced in channel height direction. Therefore a correction depending on the particle diameter-to-channel height-ratio (dP/h) is made. The bases for this correction are among other things CFD simulations, which show the temperature profile of a particle in different dP/h-ratios. The results show that only with particles showing a low dP/h-ratio an undisturbed temperature profile can be formed.

The heat flow which exists at the interface particle/melt in case of a pure thermal conduction is corrected by a thermal conduction correction factor $f_lh$.

$$f_{lh,sim} = 1,30525 + 1,98091 \cdot \frac{d}{h}$$

The correction factor is valid in the range 0,2 < d/h <0,9 and is based on an approximation of the dimensionless temperature field. It shows a good match with the CFD results, and compares the middle temperature gradient at the particle surface with finite channel height to the middle temperature gradient with infinite expansion.

An analytical solution of the energetic differential equation is in this case hard to derive [Pap06].

The factor $f_{lh}$ does not include a consideration of the expansion in channel width direction. In this case numerics are also used because an analytical derivation does not seem possible. The simulations show that the influence of the channel width only occurs with low pitch-screw diameter-ratios. But these are rarely used in practice. A correction for the heat flow with reference to the finite channel width is therefore neglected.

The influence of the convection on the melt process in twin screw extruders was until here neglected. That is why a factor fk to consider the convective heat transfer is introduced in the modified melting model. This factor is also based on CFD simulations. Based on a single rotating particle surrounded by melt in an unwound twin screw channel different relations of the pellet diameter to the channel height dp/h as well as different Péclet numbers were used. The latter balances the ratio of the convection and the thermal conduction and describes the amount of the convection occurring in the flow.

$Pe = \frac{\rho_m \cdot c_{pm} \cdot v_{0z} \cdot d_p}{\lambda_m}$

The simulation results show that for high convective shares a considerably faster melting is to be expected. Mathematically, the correction factor can be described by:

$f_k = 1 + \frac{(b_1 \cdot \kappa + b_2 \cdot \kappa^2 + b_3 \cdot \kappa^3 + b_4 \cdot \kappa^4) \cdot Pe}{(b_5 \cdot \kappa + b_6 \cdot \kappa^2 + b_7 \cdot \kappa^3 + b_8 \cdot \kappa^4) \cdot \left(\frac{3}{4} \cdot \sqrt{Pe} + b_9\right)}$

with: $\kappa = 1 - e^{-\frac{dp}{h}}$

and it compares the heat flow from the analytical calculation (without considering the convection) with the heat flow from the simulations. Therefore the following is true:

$\dot{q}_r = \left(-\lambda_m \frac{\partial T}{\partial r}\right) \cdot f_k \cdot f_{lh}$

The regression constants b1 to b9 can be taken from the following table.

Table: Regression constants

The melting process can not be considered without taking into account the whole process. Interactions between single particles as well as the particle dimension have an influence on the viscosity of the melt-solid mixture.

According to Potente and Melisch [PM96], [Mel98] the interactions are considered by using an effective channel height and width which is dependent on the amount of solid. As shown in the illustration below the melt-solid mixture is considered independent of each component. In this way the flow can be modelled through the melt above the solid layer.

Figure: Consideration of the solid particle amount in the original disperse melting model [Thü08]

In the modified melting model the flow in the melt-solid mixture is formed with the help of an adjustment of the viscosity. For this purpose the correction of the power law consistency with the factor $f_Φ$ is in introduced.

$$K_{sm} = K \cdot f_\phi$$

$$f_\phi = \frac{\eta_{sm}}{\eta_0} = 1 + \frac{5}{2} \cdot \phi_{V,s}$$

The correction factor $f_Φ$ is dependent on the volume content of the solid $Φ_{V,s}$ and compares the viscosity of the fluid-solid mixture $η_{sm}$ with the viscosity of the pure fluid $η_0$. This correction is based on Einstein but for large solid shares the viscosities calculated deviate too much. A further development according to Guth and Simah yields better results. Therefore the following is true for the correction factor:

$$f_\phi = \frac{\eta_{sm}}{\eta_0} = 1 + \frac{5}{2} \cdot \phi_{V,s} + \frac{141}{10} \cdot \phi_{V,s}$$

A conversion of the volume content on the mass content yields:

$$f_\phi = 1 + \frac{5}{2}\left(\frac{\phi_{M,s}}{\phi_{M,s} + \frac{\rho_s}{\rho_m}(1 - \phi_{M,s})}\right)\phi_{V,s} + \frac{141}{10}\left(\frac{\phi_{M,s}}{\phi_{M,s} + \frac{\rho_s}{\rho_m}(1 - \phi_{M,s})}\right)^2$$

If the $d_p/h$-ratio approaches 1 the particle diameter considered corresponds approximately to the channel height, and the influence of the particle size on the flow has to be considered. Following Pape [Pap06] another correction of the power law consistency is introduced.

$$K_{sm} = K \cdot f_\phi \cdot f_{dh}$$

$$f_{dh} = 1 + \frac{\frac{d_p}{h} \cdot \left(1 + \left(\frac{d_p}{h}\right)^2\right)}{\left(\frac{\rho_s}{\rho_m} - \frac{1 - \Phi_{M,s}}{\Phi_{M,s}}\right) \cdot \left[\left(\frac{d_p}{h}\right)^2 + \frac{d_p}{h} + 4\right] \cdot \left(1 - \frac{d_p}{h}\right)}$$

Based on a solid particle in a Newtonian melt the flow conditions around this particle are considered in order to determine the correction. After this the viscosity for the flow of an equivalent volume flow of pure melt is determined. The ratio of the viscosity of this consideration and the viscosity of the solid-melt mixture leads to the correction with the factor $f_{dh}$.

The physical-mathematical description of the modified disperse melting behavior is basically based on the assumptions of the prior melting model for disperse fillers. Here the following additional aspects are taken into consideration:

• The influence of the particle dimension on the flow $f_{dh}$,
• The influence of the particle interactions on the flow $f_\Phi$,
• The Influence of the finite expansion of the channel in height direction on the heat flow in radial direction $f_{lh}$
• The influence of the convection on the heat flow in radial direction $f_k$.

by the assumptions made.

With the influence of the correction factors flh and fk the heat flow in radial direction is given by the following equation (cf: Calculation_of_the_solid_bed_reduction_along_the_melt_path):

$\dot{q}_r = -\lambda_m \cdot \frac{\partial T}{\partial r} \cdot f_k \cdot f_{lh}$

These factors are used in the whole calculation process which was introduced in the previous chapter. Therefore the following is true for the change of the particle radius in the interval considered:

$r_{i+1} = \sqrt{r_i^2 - \frac{2 \cdot \lambda_m \cdot f_{lh} \cdot f_k}{c_m \cdot \rho_s \cdot \bar{v}} \cdot \ln\left(1 + \frac{c_m \cdot (T_m - T_{flow})}{\Delta h \cdot f_{lh} \cdot f_k}\right) \cdot \Delta z}$

The factors fdh and fΦ, which reflect the influence of the particles on the flow, are anchored in the calculation of the middle flow velocity $\bar{v}$.

Since SIGMA 11.1, a new melting model is available - the 2D disperse melting model. The naming is based on the two-dimensional particle temperature analysis - from the feeding zone through the solid-conveying-zone into the melting zone, the temperatures of the particles are calculated along the particle radius and along the screw length.

The previous melting models remain and can still be selected. However, the recommended calculation setting is the choice of the 2D disperse melting model when polymer particles with a diameter greater than 2 mm are processed.

The new model complements the modified disperse meltdown model (MDA) from 2008. One of the extensions of the MDA is the analytical calculation of the particle temperature development from the feeding zone to the melting zone. In previous calculations, only ambient air heating had an influence on the temperature of the solid particles, whereby the calculation of the air temperature represents a simple averaging. In the current 2D disperse melting model, convection from ambient air, contact with the heated barrel wall and frictional heating through friction between granulate and granulate, and granulate and steel are taken into account. For this purpose, these components are considered as heat flows in the particle. The air temperature is calculated according to VDI Wärmeatlas as heating in the annular gap and averaged over the screw length.

In various sources of literature, the heating caused by a plastic deformation in the kneading disks is the decisive influencing factor for the melting in twin-screw extruders. These findings are confirmed by suitable experimental investigations. Based on this motivation, the existing melting model was extended so that these deformations and the resulting initial partial melting could be taken into account. For this purpose, a partially analytical model has been set up, which calculates in the first kneading disks of a melting zone that particle fraction which is deformed in the intermeshing region of the twin-screw channel. Depending on the material, mass throughput, input temperature and rotational speed, this proportion undergoes a previously empirically determined energy input and thus an increase in temperature.

In order to be able to make statements about the degree of melting from the information on the temperature development in the solid-conveying and melting zones, a shell model is built up using the finite difference method and the determined heat flows in the particles: The polymer particles, which inlet temperature and diameter are known, are divided into a defined number of shells, so that for each shell a temperature depending on the initial particle temperature, the prevailing heat flows and the heating time can be calculated. Thus, over the solid-conveying and melting zone a shell temperature profile results. If one shell temperature exceeds the glass transition temperatures or crystallite melting temperatures, then the material in this shell is considered as molten material. By means of a volumetric calculation, it can thus be determined how many percent of the particle has already passed into the molten state - the degree of melting and the remaining solids content are thus recorded and transferred to the result output. Within the zone in which energy is introduced into the particle by means of deformations, the shell temperatures are increased by the temperature rise calculated within the deformation model. The addition of a temperature vector to each shell is due to the fact that no externally flowing heat flow is present, but the complete particle is heated by deformation. The procedure for determining the degree of melting can be carried out as described.

The energy input by deformation decreases, the higher the temperature of the material immediately before the deformation is. In addition, with increasing melt content and concomitant increase in the degree of filling, the heat conduction from the melt dominates the particle temperature development. It is no longer negligible. The melting of the spherical particles in the surrounding melt is well imaged by the modified disperse melting model. Therefore, this model is used when the energy input due to deformation becomes negligible.

The filling degree of the components melt and solid is the decisive criterion for going on with the calculation according to the modified disperse melting model. In this model, the spherical particles in the screw channel after the „deformation zone“ are almost closely spaced packed. It is model assumption the disperse melting starts when the voids between the spheres are filled with melt. The packing density of a densest packing is

$$\frac{\pi}{3\sqrt{2}} \sim 0,74048 = 74,048$$

which means that the free volume occupies 25.952 %. As soon as the degree of molten material exceeds the melt rate of 25.952 %, the subsequent melting is calculated by the modified disperse melting model. This modeling takes into account the melting by a successive reduction of the radius due to convective heating in a finite channel geometry.

[Ans93] Ansahl, J.: Grundlagen für die Auslegung dichtkämmender Gleichdrall-Doppelschneckenextruder, Dissertation Universität Paderborn, 1993

[Mar84] Martin, H.: Wärmeübertragung in Feststoffpartikeln, Technische Mitteilungen, (1984)12

[Mel98] Melisch, U.: „Grundlagen zur Simulation des Förder- und Plastifizierprozesses dichtkämmender Gleichdrall-Doppelschneckenextruder“; Dissertation; Universität Paderborn; 1998

[NN88] VDI - Wärmeatlas, VDI-Verlag GmbH, Düsseldorf, 1988

[Pap06] Pape, J.: „Grundlagen der Prozesssimulation von Einschneckenkonzepten zur Hochleistungsplastifizierung„; Dissertation; Universität Paderborn; 2006

[PM96] Potente, H.; Melisch, U.: „Theoretical and Experimental Investigations of the Melting of Pellets in Co-Rotating Twin-Screw Extruders“; International Polymer Processing; Volume 11; 1996; S.101-108

[Pot90] Potente, H.; Mitarbeiter: Rechnergestützte Extruderauslegung - Kunststofftechnisches Seminar an der UNI-PB, 1990

[Sch84] Schlünder, E.H.: Heat transfer to packed and stirred beds from the surface to immersed bodies, Chem. Eng. Process, 18(1984)

[Thü08] Thümen, A.: „Untersuchung und Beschreibung des dispersen Aufschmelzens in Gleichdrall-Doppelschneckenextrudern„; Dissertation; Universität Paderborn; 2008

[TM87] Tsotsas, E.; Martin, H.: Thermal Conductivity of a Packed Beds, a Review, Chem. Eng. Process, 22(1987)