1. Introduction
Single-screw extrusion is one of the most widely used production processes in the plastics industry. Semi-finished products such as tubes, hoses, foils, plates, and geometrically even more complex profiles can be produced [1,2]. Hence, single-screw extrusion is still the subject of current industrial and fundamental research. In recent years, computer-aided modeling of single-screw extrusion and research for new simulation approaches has become increasingly important. In general, the existing studies can be subdivided into investigations of the feed and the melting section. First attempts to describe mathematically the transport of solids in the feed section in one dimension go back to Darnell and Moll [3] and Schneider [4], who made the basic assumption that the interparticular friction is greater than the friction between particles and walls. Schneider [4] formulated the first equations for calculating flow rate and pressure build-up. Over the years, the first approaches were further developed in [5,6,7], until Schneider [8,9] suggested the insertion of axial grooves in the feed section in order to convey independently of pressure. Further studies [10,11,12] concentrated on the modeling and calculation of grooved feed sections.
In addition to Potente and Pohl [13], Moysey and Thompson [14] were pioneers in 2004, when they applied for the first time a three-dimensional numerical approach, the so-called discrete element method (DEM), to a smooth feed section. This method has the advantage that each granule and its interaction with surrounding particles and walls can be mapped individually. However, potentially higher computing times are disadvantageous. Over time, the results obtained by Moysey and Thompson [14,15] and the DEM itself proved to be very promising for calculating the feed section. Since then, further works such as Schöppner et al. [16,17] and Celik and Bonten [18] focused on the application of DEM. In addition, Leßmann [19,20] investigated a smooth extruder with rotational speeds up to 900 min−1 and considered axial forces acting on granules.
The DEM also proved suitable for modeling the transport of solids in twin-screw extruders. Interesting applications can be found in Carrot et al. [21] and Potente et al. [22]. Besides the solids transport itself, Schöppner [23] also modeled the heating of granules along a twin-screw extruder.
However, despite the valuable work done in recent years, the models are partly dependent on empirical model parameters. What is much more important though is that the investigation by means of DEM is always limited to the feed section, since the phase transition during melting cannot be modeled. As a result, the melting section, limited by the available numerical simulation possibilities, is considered largely separate from the feed section.
Initial investigations of the melting behavior go back to Maddock [24], who stopped an extruder during operation, let it cool down and then pulled the screw out of a barrel. He found that a melt film forms between a granular bed and a cylinder wall. Based on these investigations, Tadmor [25,26] developed a mathematical model to describe the melting behavior. Donovan [27] introduced extensions in terms of heating a granular bed along the screw, whereas Pearson [28] later assumed a nonconstant melt film thickness. Potente [29] also introduced analytical equations based on the investigations of Tadmor [25,26] to calculate the melting behavior.
Besides the analytical–empirical approaches mentioned so far, Viriyayuthakorn and Kassahun [30] pursued the aim to describe the melting process by means of conservation equations and equations of state, in order to investigate a model that is not limited to materials and process parameters. The melting process is not described by a phase change, but rather by the associated change of material properties. They chose the specific heat capacity as a decisive parameter. However, an experimental validation was not carried out. Altinkaynak’s [31] approach, based on the finite element method (FEM), was also basically developed out of this idea. However, Altinkaynak [31] chose the viscosity that changes during melting as a parameter. For the solid state, the viscosity is assumed to be infinitely high and decreases with increasing temperature. Hopman et al. [32] chose an approach based on the finite volume method (FVM) with a similar idea.
A very good overview of the basics of global modeling of single-screw extrusion and the existing modeling approaches can be found in Wilczyński et al. [33]. The authors concluded that a coupled approach via Computational Fluid Dynamics (CFD) and DEM for the holistic modeling of single-screw extrusion proves promising.
However, it is also required that a suitable melting model can precisely describe the phase transition. Karrenberg et al. [34] also criticized about the lack of a modeling approach combining CFD and DEM.
The present work resolves this aspect. The current state of research leads to the conclusion that the holistic modeling of single-screw extrusion is only possible under separate consideration of the respective feed and melting section. In the authors’ view, this is because until a few years ago there was simply no simulation tool that allowed describing the interaction of particles and fluids in three dimensions. Modeling the coexistence of several phases remains a challenge to this day, which is reflected in the fact that solid and liquid states are considered decoupled from each other. Furthermore, from the authors’ point of view, a suitable melting model that describes the phase transition of plastic granules is an important aspect but is missing, thus enabling the holistic, three-dimensional simulation of the feed and melting sections to be necessary.
Therefore, the primary goal of this work is to explore a novel melting model for plastics to describe phase transitions. This melting model is to be implemented into the simulation software CFDEMCoupling® [35]. This is an open-source coupled CFD-DEM framework which extends OpenFOAM® [36] to include a coupling with the DEM software LIGGGHTS® [37], which stands for LAMMPS improved for general granular and granular heat transfer simulations.
In order to ensure reliable results, the aim is also to validate the melting model experimentally. For this purpose, a melting apparatus will be set up and put into operation. This apparatus reproduces similar physics appearing in a real single-screw extruder. With the help of the experimentally obtained data, conclusions can be drawn about the reliability of the melting model. The overriding and long-term goal is to apply the newly researched simulation method to single-screw extrusion in the sense of holistic simulations. 2. Novel Modeling Approach
The new melting model stands out from the current state of the art, as the melting behavior is not represented by a temperature-dependent change in the dynamic viscosity or the specific heat capacity (see [30,31]). The novelty of the melting model is given by a general CFD-DEM approach [38], which couples the so-called volume of fluid method (VOF) [39,40] with the DEM. The cutting-edge features of this model include the physical interaction between three phases (air, solid, and molten particles), the heterogeneous temperature distribution within a particle and the particle melting representation. Regarding the latter feature, a decrease of the radius of a DEM particle as a function of the temperature during the melting process is explicitly represented. The main requirement of the model is to describe the phase change of plastic granules under consideration of heat conduction, convection, particle heating due to friction, and the resulting mass change.
In the following, the focus lies in the description of the modeling approach. First, the existing basic equations of the DEM are introduced. Then, all extended model equations including the novel melting model are presented. 2.1. DEM Approach
For the description of particle–particle and particle–wall interactions, the basic equations of the DEM were used, which follow a Lagrangian approach. The simulation software LIGGGHTS® (version 4.1.0) was used for this purpose. In the DEM, the trajectories of all individual discrete particles were followed, and their velocity and acceleration were determined by deriving the trajectory curve according to time. According to Newton’s equation of motion, the product of massmiand accelerationx¨iof a particleiis equal to the sum of all acting forcesF⇀i(Equation (1)). The force and moment balances are shown in Equations (1) and (2), respectively:
mi x¨i=F⇀i,n+F⇀i,t+F⇀i,f
I⇀idωidt=R⇀i×F⇀i,t+T⇀i,r
whereF⇀i,nis the particle–particle contact force in the normal direction,F⇀i,tis the particle–particle contact force in the tangential direction,F⇀i,fis the resistance force, which acts on the particles through the surrounding fluid,I⇀iis the moment of inertia andωiis the angular velocity,T⇀i,r describes an additional torque that can be used to simulate nonspherical particles [38], andR⇀iis the position vector of the contact point of two particles. In addition to the forces mentioned here, other forces, such as electrostatic or magnetic forces, can also be taken into account. Furthermore, the interparticular interaction is described by means of contact models. To calculate the contact forces, the linear spring-damper model is often useful. The normal force is thus composed of the spring and damper forces:
F→i,n=−kn·δn+cn·Δu→n
whereknandcnrepresent the spring and damper coefficients,Δu→ndescribes the relative velocity of the particle in the normal direction, andδndescribes the virtual overlapping of two particles. In a real collision of two particles, a part of the energy is converted into heat by friction and plastic deformation. This dissipation energy is represented by damping elements with the damping coefficients in the normal directioncn in the model. According to Hertz and Mindlin [41,42,43],knandcncan be calculated from the material properties in the normal (indexn) and tangential directions (indext). The tangential force is limited by the static friction force and consists of the shearing force and the damping force in the tangential direction:
F→t=min{|kt·∫tctΔu→t·dt+ct·Δu→t|, |μ·F→n|}
In Equation (4),tcis the time when particles come in contact,Δu→tis the tangential component of the relative velocity, andμ stands for the coefficient of friction. In addition to the balance of forces and moments, the energy equation for discrete particles must also be considered for modeling the melting behavior. The basic equations for the description of heat conduction in LIGGGHTS® are originally based on the investigations of Chaudhuri [44]. Accordingly, the heat fluxq˙i,jcondbetween the individual particlesiandjdue to heat conduction is calculated as follows:
q˙i,jcond=hi,j·ΔTi,j
whereTi,jcorresponds to the particle temperature andhi,jto the interparticular heat transfer coefficient, which is calculated according to Equation (6) from the thermal conductivity coefficientki,jand the contact areaAconof the colliding particles:
hi,j= 4ki kj ki+kj·(Ai,j)1/2
whereki,jis surface-independent and has therefore the dimensionJK·s·m, whereashi,jhas the dimensionJK·s. The temporal temperature development for a discrete particledTidtcan be calculated according to the following energy balance:
mi cidTidt= ∑i−jq˙i,j+ q˙src/sink
wherecicorresponds to the specific heat capacity andmiis the mass of the particle. The first term on the right side of Equation (7) represents the heat flux due to particle interaction, and the second term describes heat sources and sinks. The conventional DEM is based on ideal spherical particles, which is a valid assumption for most applications.
In a real melting process, not only particle–particle and outer particle–wall heat transfer is of special interest. The rate of heat transferq˙due to convective heat transfer through the surrounding fluid to the particle must also be considered and can be calculated according to Equation (8):
q˙=hconv·A·(T∞−Ti)
whereT∞is the temperature of the surrounding fluid andTiis the temperature of the particle;hconv represents the heat transfer coefficient and is not a material parameter. It is generally difficult to determine for bulk materials and is calculated in LIGGGHTS® using the Nusselt number. More about the convective heat transfer model used can be found in Li and Mason [45].
2.2. Extended CFD-DEM Approach
For the CFD-DEM simulation, the software CFDEMCoupling® was used, which is an extension of OpenFOAM® (version 5.x), enabling coupled simulations with LIGGGHTS®. Before the advanced CFD-DEM conservation equations are presented, a schematic 2D consideration of the coexistence of three phases in any finite volume element (CFD cell) is necessary (Figure 1).
The three phases can be divided into the solid state (plastic granules here represented by discrete particles), the liquid state, which consists of the melt that forms during the melting process, and the gas phase represented by the surrounding air.
The so-called voidfractionεat a CFD cell is defined as the ratio between the volume occupied by the fluidVf(gas phase plus liquid phase) and the total volume of the CFD cellVcell:
ε=VfVcell=1−εs
whereεsis the solid fraction, defined as the ratio between the volume occupied by the solid fractionVsand the volume of the CFD cell. Hence,εassumes to be 1 if there are no particles in the CFD cell and strives towards the value of 0 if more particles are in the CFD cell.
The total fluid volumeVfresults from the sum of the volumes of the individual phases:
Vf=Vl+Vg
These must be calculated partially. The respective phase fractionαiresults from the ratio of the volume of the respective phaseVito the total fluid volumeVf:
αl=VlVf
αg=VgVf
The following volume balance can be drawn up:
Vs+Vf·(αl+αg)=Vcell
By diving both terms of the equation byVcell, the following equation is obtained:
εs+ε·(αl+αg)=1
In light of this consideration, the conservation equation for phase fractions can now be established in the form of a transport equation. The transport of the phase fractionsαlandαgtakes place via the flow velocityu⇀f. In addition, the temporal mass changes˙meltmust be considered due to the phase transition during melting. The mass conservation equations are then shown as follows:
∂(ρl αl)∂t+∇·(ρl αl u⇀f)=s˙melt
∂(ρg αg)∂t+∇·(ρg αg u⇀f)=0
Since the mass of the gas phase is preserved and is not converted, the right side of Equation (16) is zero. The conservation of momentum (see Equation (17)) is analogous to an unresolved CFD-DEM approach according to Kloss et al. [38]. With an unresolved CFD-DEM approach, it must be ensured that each CFD cell in the computational area provides at least the volume of one spherical particle.
∂(ρf αf u⇀f)∂t+∇⋅(ρf αf u⇀f u⇀f)=∇⋅(αfτ_)−αf∇p+αf ρfg⇀+σκ∇(αl)|∇(αl)|−Ksf(u⇀f−〈u⇀p〉)
In Equation (17),αfis the fluid volume fraction occupied by the liquid–gas phase,∇pis the pressure gradient,τ_is the stress tensor,αf ρfg⇀is the gravitational term,σindicates the surface tension between the gas and liquid states,κis the curvature of the gas–liquid interphase, and the expressionKsfdescribes the exchange of moments between the fluid and solid phases taking into account the particle velocityu⇀p . In the past, there have been various approaches to modeling the momentum exchange term [46,47,48]. A combination according to Gidaspow et al. [49] of the Ergün equation [50] for fluid volume fractionsεf of ≤0.8 and a model according to Wen and Yu [51] forεf>0.8 has become generally accepted. For further details, please refer to Kloss et al. [38].
The transport of the fluid phase volume fractionαf is based on the VOF method [39,40] and is described by Equation (18):
∂(αl)∂t+∇⋅(αl αg u⇀f)−∇⋅(αl(1−αl)⋅u⇀c)=−αl∂(αg)∂t
where the source term on the right side accounts for the fluid displaced by the particles. According to Rusche [52], the compression term including the compression velocityu⇀cleads to a sharper interphase between the two different fluid phases.
2.3. Novel Melting Model for Spherical Solids In order to describe a phase transition from the solid state to melt and the resulting change of the mass over time, a novel melting model was developed and implemented into the simulation environment CFDEMCoupling®. This melting model mathematically describes the unsteady heating of a discrete spherical particle using the unsteady heat conduction equation. In reality, the heating of a plastic granule takes place from outside to inside. This results in an inhomogeneous temperature field. In the model presentation, a spherical particle is therefore discretized for the mathematical description of the temperature distribution into “shells”. One spherical particle is therefore composed of any number of shells.
As an example, in Figure 2a, it is assumed that a particle is composed of three shells. If a particle is heated until the temperature of the outermost shell (T3) reaches the specified melting temperature (Tmelt), the outermost shell enters the melt phase. Two more shells then remain, which also enter the melt phase as soon as their temperatures (T2,T1) exceed the melting temperature.
It should also be noted that the melting of the shells is accompanied by a reduction of the particle radius. As soon as the particle radius falls below a specified minimum valuermin, the particle is considered completely molten.
The unsteady temperature distribution within a shell is calculated according to Equation (19) using the one-dimensional unsteady heat conduction equation in spherical coordinates, assuming that only a temperature gradient is formed in the radial directionr:
ρs cp∂T∂t=λs(2r∂T∂r+∂2T∂r2)
For the calculation of the temperature distribution, the thermal conductivity of the particleλs, its densityρs, and specific heat capacitycpmust therefore be known. In a good approximation, these material properties can be regarded as temperature-independent in the respective temperature range. To calculate the temperature distribution within a shellkforN shells (Figure 2b), Equation (19) is linearized using a backward Euler method:
Tn+1, k−Tn, kΔt=D(Tn+1, k+1−Tn+1, k−1rΔr+Tn+1, k+1−2Tn+1, k+Tn+1, k−1Δr2)
wherencorresponds to the discretized time step, the termλsρs cpwas combined to form the diffusion coefficientD. For the solution of the one-dimensional unsteady heat conduction equation, a boundary condition for the calculation of the heat fluxqwith the environment at the positionr=Ris furthermore given as:
q|r=R=− λs∇T
With the equations listed so far, the temperature distribution in a particle can be calculated that exchanges heat with other particles in the environment with a wall or with a surrounding fluid. For the melting process, the change of the particle radius and especially the time-varying molten mass must also be taken into account.
The melting rate is calculated according to Equation (22) and is transferred for the CFD calculation to Equation (15):
s˙melt=mmeltΔt=−(−λs·C+β)L ·A
whereLis the mass related melting enthalpy, which is determined experimentally by means of differential scanning calorimetry (DCS),βdescribes the heat flux through the particle surface. In the absence of sources or sinks, this term corresponds to the termq|r=Rin Equation (21), related to the surface of the sphere,Ais the surface area of the discrete spherical particle through which the heat exchange occurs, andCis equal to the temperature gradient between two shells and described as:
C=Tn, k−Tn−1, kΔr
The reduction of the particle radius due to melting is described by the termΔrsimilar to Equation (24):
Δr=RN−k
whereRis the initial radius of the particle andN−kdescribes the decreasing number of shells within a particle depending on time and temperature.
2.4. Calculation of Frictional Heating
For smooth and spherical particles, the dissipation and frictional heatPdissin the normalnand tangentialtdirections can be calculated via the respective damping forceFand the relative velocity componentun,tat the contact point:
Pdiss=Fn,t·un,t
In case of contact between two particles, the power dissipated is distributed among the two as follows [53]:
Pdiss, 1=12·(k1k1+k2)·Pdiss
Pdiss, 2=12·(k2k1+k2)·Pdiss
wherekrepresents the thermal conductivity of the generic particle. The model employs Equations (25)–(27) to calculate the respective heat flux for each individual particle.
3. Melting Apparatus and Simulation Model 3.1. Principles and Functionality
In order to validate the melting model presented in Section 2.3, a novel melting apparatus was set up and put into operation within the scope of this work. The apparatus is based on the work of Chung 2010 [54] and is intended to determine the melting performance of a real single-screw extruder by scaling up the screw diameter. The extruder conditions that contribute to the melting process are imitated in the melting apparatus. Similar devices are called “screw simulator” in the literature [54].
The schematic sketch in Figure 3 illustrates how it works. A stamp with the cross-sectional dimensions of 50 mm × 30 mm, driven by a hydraulic system, compresses in a prechamber by its axial movement the previously inserted and weighed quantity of plastic granule. The force of the stamp is adjustable, and the axial position is detected by a distance-measuring system. The maximum length covered by the stamp is 300 mm.
The prechamber can be heated up to 200 °C, but a significantly lower temperature is set during the experiments to prevent premature melting of the granules. The actual melting is performed by the rotating ring which is heated to 200 °C by means of a slip ring and is located directly downstream of the prechamber. The ring’s shape can be designed as desired to investigate geometric factors influencing the melting behavior. For example, a smooth surface or an axially or helically grooved surface can be used. This means that the rotating ring takes over the task of the barrel in single-screw extrusion. The ring has a diameter of 300 mm and is mounted on a rotating shaft driven by an electric motor. The shaft is supported against axial displacement by spherical roller bearings. The surrounding housing is heated to ensure almost adiabatic conditions. The melt film that forms on the surface of the ring can be removed with a melt scraper. The mass of plastic granules molten over time is then weighed.
Thus, by means of the melting apparatus, the melting rate can be determined experimentally under different process conditions such as speed and temperature for different materials. For the digital signal processing and output of the data a user interface within the used software LabVIEW® [55] was developed.
3.2. Process Window
The melting rate was the target value of the experimental investigation, as it served as a comparative value for the evaluation of the melting model. As a parameter, the rotational speed was varied according to the overview in Table 1. All the other variables such as the stamp the force with which the granules were compacted and the temperatures in the housing in a prechamber and in a ring were kept constant during the test. The rotating ring was grooved. The 17 grooves had a width of 5.5 mm. For each operating point, the melting rate was determined six times to ensure sufficient accuracy. The procedure for the experimental determination of the melting rate can be described as follows:
- Predrying of the plastic;
- Heating up the melting apparatus;
- Filling the prechamber with a defined amount of granules when the stamp is moved out;
- Control of the stamp and start of the drive which sets the ring in a rotational movement; and
- Removal of the forming melt film and stopping the time as soon as a steady state is reached.
3.3. Materials and Characterization
A polypropylene, Moplen RP 210G [56] by LyondellBasell, Rotterdam, Netherlands, was used for the tests. The mass of the granules filled into the prechamber was always 100 g. In the following, the material used was referred to as PP.
For the performance of CFD-DEM simulations, a large number of material parameters were required. Since the melting process was simulated in this work, thermodynamic parameters were required in addition to the mechanical parameters that are usually necessary for a pure DEM simulation. To achieve reliable simulation results, these should be determined as accurately as possible and included into the simulation model.
Table A1 (Appendix A) gives an overview of some important material parameters and the respective determination method or equipment used. For example, to determine Young’s modulus, tensile bars according to ISO 527-2:2012—Type 1A were injection molded on an Arburg Allrounder 320 S from the manufacturer Arburg, Lossburg, Germany. Furthermore, the particle size distribution was determined with a modern scanner and was also considered in the simulations.
3.4. CFD-DEM Modeling of the Melting Apparatus
The melting apparatus was modeled according to the experimental setup. In general, a CFD-DEM simulation should be limited to the respective domains. With the apparatus shown in Figure 3, the melting process is expected to take place in the prechamber. The housing surrounding the rotating ring plays a subordinate role, since the largest part of the melt that forms is discharged via a melt scraper. For this reason, it is sufficient to limit the modeling and simulation of the apparatus to the prechamber, the stamp, and the grooved ring (see Figure 4). This also saves computing time, because the computing area is considerably reduced.
For the calculation, a basic distinction must be made between the types of computational meshes. In CFD simulations volume meshes were required, whereas surface meshes were required in DEM simulations. Both grid types were generated with the help of ANSYS ICEM CFDTM [57].
On the DEM side, the compaction of a granular bed was modeled by the axial movement of the stamp over a simple rectangular surface. The cross-section corresponded to the dimensions of the stamp of the experimental setup. Furthermore, the rotation of the grooved ring and the resulting force transmission to the particles were taken into account. The grooves were completely resolved. The stamp and the grooved ring were imported into LIGGGHTS® as surface meshes in the Standard Triangulation Language (STL).
The CFD domain represents the prechamber and the rotation of the grooved ring. This is shown in blue in Figure 4. As there is no flow, it only represents an air-filled volume in this case. However, this volume is essential, under which the calculation of the melting rate is possible according to the presented CFD-DEM equations.
The temperature boundary conditions were specified both in the CFD and in the DEM according to the experiments (see Table 1).
As explained in Section 2.2, the CFD computational mesh generation in an unresolved CFD-DEM approach represents a special challenge, since there is a dependency between the volume of the CFD cell and the particle volume. The cell size is thus limited by the particle diameter. Therefore, a computational mesh study was performed in advance to ensure mesh and time convergence. More detailed explanations including the results can be found in Figure A1 (Appendix B). The software ParaVIEW [58] was used to visualize the simulation results.
4. Results and Discussion
The discussion of the results focuses on the evaluation of the melting rate from the experiment and the simulation and their comparison. First, in Section 4.1, the influence of the number of shells on the simulation result is highlighted. Afterwards, in Section 4.2, the melting rates obtained using the new melting model are compared and discussed with experimental results of the measuring apparatus.
4.1. Influence of the Number of Shells
After the successful implementation of the melting model into the CFDEMCoupling® framework, the first issue to be clarified is how the number of shells per DEM particle influences the melting behavior. For this purpose, CFD-DEM simulations using PP were performed for the melting apparatus and for a number of shells (see Figure 2) ranging from 1 to 100. In order to save computing time, the rotational speed was set to zero. The stamping force was 100 N. The percentage of the molten massξcan then be determined by evaluating the simulations as follows:
ξ=(mtotalm−1)·100
wheremtotalis the total mass specified in the simulation at time 0 and it was deliberately set to 12 g in the simulation to reduce computational time, which was reasonable, since the change of mass per time was of interest;mstands for the value of the mass at any further time and was returned from the simulation.
Figure 5 shows the course ofξfor different numbers of shells over a simulation time of 60 s. A shell number of one corresponds to the case that the whole particle was considered as molten when the particle temperature exceeded the melting temperature, i.e., no discretization of the particle took place.
All the characteristic curves are similar: Up to approximate 10 seconds, the particles were heat up andξwas approximately 0. In the further course,ξrose continuously. It can be observed that by increasing the number of shells from one to three, the start of melting was shifted. This can be explained as follows: It took less time to melt the outermost shell of a particle with three shells than a complete, nondiscretized particle (shell number = 1). The outermost shell of a particle with three shells thus melted earlier, which was reflected in the increase in the molten mass.
This means that three shells melted up more in percentage terms and the gradientdξdtwas slightly steeper. If the number of shells was further increased to 10, 50, and 100, no significant influence can be seen, demonstrating that the discretization achieved sufficient accuracy with three shells.
4.2. Comparison of the Calculated and Measured Melting Rates
First, the CFD-DEM simulation results based on the implemented melting model are illustrated. The left column in Figure 6 shows the melting phenomena in the DEM domain of the melting apparatus. The stamp compresses the particle bed in the prechamber with a force of 1000 N and presses the particles onto the grooved ring, which rotated at a speed of 3 min−1. In the left column, the reductions of the particle radius and particle number due to phase transition are shown. The right column in Figure 6 shows the melting phenomena in the prechamber and the grooved ring in the CFD domain.
For reasons of clarity, a plane was defined exactly in the middle through the prechamber. At the same time in the process as the DEM simulation, the formation of melt in the prechamber can be observed on the basis of the corresponding phase fraction.
At time zero, the particles were present in their natural size distribution in a loose bulk. At this time, no melt formed in the CFD domain, as can be seen in the right column of Figure 6. After 16 s, the particles were already compressed and partially molten. This can be seen from the fact that the volume fraction of the melt immediately above the grooved ring was greater than zero. The melting process therefore took place as expected on the heated grooved ring. This is where the highest temperature occurred. The melting process was also greatly enhanced by the frictional power of the grooved ring.
In the further course of the coupled simulation, the number of particles decreased steadily, and melt continued to form, which was afterwards transported. Similar behavior was observed for the other considered operating points at rotational speeds of 2 and 3 min−1.
After the melting process was illustrated, a comparison between the results of simulation and experiment was drawn. Figure 7 shows the measured melting rates for three rotational speeds compared to the simulation results.
The melting rate increased with increasing speed under otherwise constant conditions. This observation is in line with the existing state-of-the-art investigations [31,59] and thus confirmed the correct functioning of the melting apparatus. It can be seen that the simulated melting rates were lower than the experimental values for all the rotational speeds.
The deviations of simulation results from the experiment were between 1.68% at 1 min−1 and 21.71% at 3 min−1. The deviations from the experiment can be attributed to the following unavoidable simplifications of the melting model and the limits of the DEM:
- The melting model was set up for ideal spherical particles. However, the PP used in the experiment had an ellipsoidal granular form.
- In the experiment, granules were deformed by the stamp force and the rotation of the grooved ring, favoring the melting process. However, the DEM cannot completely map this kind of deformation of the particles. That is why the simulation values were lower than the experimental results.
- The Young’s modulus, which has a considerable influence on the melting behavior, was in reality temperature-dependent, but was assumed to be constant in the DEM. A temperature-dependent implementation of the Young’s modulus is not possible at present, as the Rayleigh time and the numerical stability of the simulation are directly influenced.
The deviations may seem large at first glance. However, the results achieved with the new modeling approach are very promising and overall in good agreement with experimental results. In regard of the evaluation of the results, it is even more important to mention that no parameter or model adjustments were undertaken in the coupled simulation of the melting apparatus. Only experimentally determined substance data were used instead. 5. Conclusions and Outlook Within the scope of this study, a novel melting model was developed. This describes the phase transition from the solid state of the plastic to the melt state by means of the CFD-DEM approach for the first time. The melting model considers the temperature change due to heat conduction, convection, and frictional heating. For the validation of the simulation results, a melting apparatus was set up and put into operation. Melting experiments were performed to determine the melting rate of PP. In parallel, the experimental setup was modeled. The influence of the number of shells on the particles’ melting behavior was drawn. Considering the necessary simplifying assumptions in the model and limitations in the CFD-DEM approach, a comparison between the simulation and the experiment showed a very good agreement without any adjustment of the parameters. To confirm the transferability of the melting model, the experiments were extended to further experimental parameters and materials. In particular, the stamp force, the ring temperature, and the geometry of the ring should be varied. As a next step, the model should be transferred to the single-screw extrusion. The modeling of the extrusion is currently proving to be much more challenging, since dynamic computational meshes must be used on the CFD side in order to correctly represent the conveying process of granules. However, a CFD-DEM simulation model of the extrusion has been successfully created.
Parameter | Value |
---|---|
rotational speed | 1, 2, and 3 min−1 |
housing temperature | 200 °C |
ring temperature | 200 °C |
prechamber temperature | 130 °C |
stamp force | 1000 N |
rotating ring | helically grooved |
Author Contributions
Conceptualization, A.C.; methodology, A.C.; software, C.K., C.G. and R.T.; validation, A.C.; simulation results, A.C.; experimental investigation, A.C.; writing of the original draft preparation, A.C.; writing of review and editing, A.C., C.B., C.K., C.G. and R.T., visualization: A.C., supervision, C.B. All authors have read and agreed to the published version of the manuscript.
Funding
This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within the DACH Lead Agency agreement (grant number: 1600/31-1, project number: 324934383).
Institutional Review Board Statement
Not applicable.
Informed Consent Statement
Not applicable.
Data Availability Statement
Not applicable.
Acknowledgments
Special thanks go to the company LyondellBasell Industries, who provided the materials used for experiments free of charge.
Conflicts of Interest
The authors declare no conflict of interest.
Appendix A. Overview of Important Required Material Data
Table
Table A1.Overview of the required parameters for the CFD-DEM simulation and determination method.
Table A1.Overview of the required parameters for the CFD-DEM simulation and determination method.
Mechanical and Other
Properties | Method of Determination and Equipment Used |
---|---|
Young's modulus | Tensile test acc. to DIN EN ISO 527 (temperature-dependent) |
Coefficient of friction | Literature value |
Poisson's ratio | Literature value |
Coefficient of restitution | Universalprüfmaschinen ZPM1455, Zwick GmbH, acc. to ISO DIN EN ISO 527 |
Density (solid state) | Data sheet PPRP210G |
particle size distribution | High-resolution scanner Epson V850 pro, Epson K.K. |
Thermodynamic properties | |
Specific heat capacity | Differential scannig calorimetry acc. to DIN EN ISO 11357-4 |
Thermal diffusivity | Laser flash analysis |
Latent heat | Differential scannig calorimetry acc. to DIN EN ISO 11357-4 |
Density (liquid state) | Acc. to DIN EN ISO 1183-1 |
Appendix B. Mesh Convergence Study
With a coupled CFD-DEM simulation, the sole analysis of the mesh influence and the given CFD time step on the simulation results was not as sufficient as with a pure CFD simulation. In addition, the influence of the selected DEM time step must be considered, since the solver does not adapt the time step automatically. For this investigation, the melting apparatus was simulated.
To save computing time, the simulation time was set to twelve seconds, whereby the initialized mass of the particles was 0.012 kg. The molten mass was evaluated for two different computational meshes and DEM time steps (see Figure A1).
Polymers 13 00227 g0a1 550
Figure A1.Molten particle mass depending on the number of calculation elements and the selected DEM time step.
Figure A1.Molten particle mass depending on the number of calculation elements and the selected DEM time step.
Figure A1 shows that the mesh and time step convergence was given for the coarser computational mesh with 7.232 elements at a DEM time step of 5 × 10-6 s. The CFD time step was chosen two orders of magnitude higher than the DEM time step and was thus 5 × 10-4 s. The average relative deviation due to the grid influence was less than 1%. Within different time steps, the average relative deviation was 2.17%. It was therefore sufficiently small to confirm the time-step convergence.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2021. This work is licensed under http://creativecommons.org/licenses/by/3.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Abstract
Existing three-dimensional modeling approaches to single-screw extrusion can be classified according to the process sections. The discrete element method (DEM) allows describing solids transport in the feed section. The melt flow in the melt section can be calculated by means of computational fluid dynamics (CFD). However, the current state of the art only allows a separate consideration of the respective sections. A joint examination of the process sections still remains challenging. In this study, a novel modeling approach is presented, allowing a joint consideration of solids and melt transport and, beyond that, the formation of melt. For this purpose, the phase transition from the solid to liquid states is modeled for the first time within the framework CFDEMCoupling®, combining CFD and DEM by a novel melting model implemented in this study. In addition, a melting apparatus for the validation of the novel melting model is set up and put into operation. CFD-DEM simulations are carried out in order to calculate the melting rate and are compared to experimental results. A good agreement between the simulation and experimental results is found. From the findings, it can be assumed that the CFD-DEM simulation of single-screw extruder with a joint consideration of the feed and melt section is feasible.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer