Discrete Element Method Modeling of the Rheological Properties of Coke/Pitch Mixtures

Rheological properties of pitch and pitch/coke mixtures at temperatures around 150 °C are of great interest for the carbon anode manufacturing process in the aluminum industry. In the present work, a cohesive viscoelastic contact model based on Burger’s model is developed using the discrete element method (DEM) on the YADE, the open-source DEM software. A dynamic shear rheometer (DSR) is used to measure the viscoelastic properties of pitch at 150 °C. The experimental data obtained is then used to estimate the Burger’s model parameters and calibrate the DEM model. The DSR tests were then simulated by a three-dimensional model. Very good agreement was observed between the experimental data and simulation results. Coke aggregates were modeled by overlapping spheres in the DEM model. Coke/pitch mixtures were numerically created by adding 5, 10, 20, and 30 percent of coke aggregates of the size range of 0.297–0.595 mm (−30 + 50 mesh) to pitch. Adding up to 30% of coke aggregates to pitch can increase its complex shear modulus at 60 Hz from 273 Pa to 1557 Pa. Results also showed that adding coke particles increases both storage and loss moduli, while it does not have a meaningful effect on the phase angle of pitch.


Introduction
Carbon anodes for the aluminum smelting process are made using a well-designed recipe of pitch binder and calcined coke particles. Anodes are consumed during the electrolysis process in Hall-Héroult cells, and it is estimated that around 17% of the cost of the aluminum smelting process comes from carbon anodes [1]. Thus, the properties and quality of the anodes have a direct effect on the performance and economy of the Hall-Héroult process. Anodes are made by mixing granulated coke aggregates with coal tar or petroleum pitch at 150˝C and the mixture, known as anode paste, is then compacted and baked. It is essential to have a good understanding of the viscosity and rheological properties of pitch and pitch/coke mixtures to have better control of final compacted anode properties. Micromechanical models can provide considerable information about the flow and compaction properties of anode paste.
Discrete element method (DEM) is used to simulate the behavior of granular materials in industrial applications especially where the dynamics and flow of a particulate material are of interest.
This method was introduced for the first time by Cundall and Strack [2] in 1979. In DEM simulations, rigid discrete elements, which are spheres in 3D and discs in 2D models, are used to model the granular material. The contact law between the elements defines the mechanical behavior of the bulk material.
The discrete element approach has attracted the interest of researchers in mining, civil engineering, pharmaceutical industries, and materials engineering to simulate the flow [3,4], compaction [5,6], and mechanical properties [7,8] of single-and multi-phase materials. Recent advancements in the performance and power of computers have provided new insights to simulate the mechanical and rheological properties of asphalt concretes and mastics using two-and three-dimensional DEM models. For instance, Dondi et al. [9] studied the effects of aggregate size and shape on the performance of asphalt mixes. Khattak et al. [10] used an imaging technique coupled with DEM to create a micro-mechanical model of hot mix asphalts (HMA). They reported that the DEM model provided a good prediction of the dynamic modulus and strength of HMAs. Since materials properties in DEM are modeled by assigning appropriate contact models to the elements, different elastic, elasto-plastic, and viscoelastic models are developed. For asphalt mixes, usually Burger's viscoelastic model is used for the mastic, and a simple elastic model is used for the aggregates [8,11]. Ma et al. [11] used Burger's model in a three-dimensional DEM model to investigate the effects of air voids on creep behavior of asphalt mixtures.
Burger's model can be also embedded in traditional linear elastic contact models in which normal and shear stiffness of the contact changes by time to include viscous deformations. Vignali et al. [12] and Dondi et al. [13] adapted this approach of viscoelastic modeling to predict the rheological behavior of bituminous binders.
In the present work, three-dimensional imaging is used to capture the irregular shapes of coke aggregates. Then, a 3D DEM model of pitch is developed to predict the rheological properties of pitch and coke/pitch mixtures at 150˝C.

Theory
A three-dimensional DEM model is composed of a combination of discrete spheres and walls. At the beginning, the position of all elements and walls are known so that the active contacts are easily determined. Then, according to the mechanical behavior of the material, an appropriate force-displacement law is applied to each contact, and the contact forces are calculated. Newton's second law of motion is then used to update the position and velocity of each element. One common contact model widely used in DEM simulations is the linear contact model. This model is simply defined by assigning normal and shear stiffness values to the contacting elements. Contact forces can be calculated according to the extent of the overlap (for normal contact force) and tangent movement (for shear contact force): Irregularly shaped particles can be generated as a clump composed of several touching or overlapping balls. Contact force calculations for balls within a clump is skipped during calculation cycle and only the contacts of a clump with neighboring clump/balls or walls are considered.
Burger's four-element model is the most common model used to simulate the viscoelastic properties of mastics and binders [14]. Burger's model, as shown in Figure 1, is composed of the Maxwell model in series with the Kelvin model; thus, it is capable of representing the material behavior under both creep and relaxation [14]. Deformation of Burger's model is the sum of the deformations of Kelvin's elements and Maxwell's elements. This means that the total deformation can be written as: Deformation of the Maxwell element, in turn, comes from the dashpot and the spring: The first and second derivatives of Equation (4) can be written as: The force-displacement equation of Burger's model can be written as: Although Burger's model is widely used by some researchers, mostly to simulate the mechanical behavior of asphalt mastic, it has not been publicly implemented in YADE [15]. Therefore, authors use an in-house implementation of this model in YADE, creating a cohesive viscoelastic model to simulate pitch. A new type of material, the so-called CohBurgersMat is defined in YADE. This material takes four of Burger's parameters in normal direction and four of Burger's parameters in shear direction. Interactions of elements having CohBurgersMat are governed by Equation (7).
Viscoelastic properties of a material are normally presented by two parameters of a complex shear modulus, G*, and phase angle, δ [13]. These parameters can be measured using a parallel-plate measuring system as schematically shown in Figure 2. The dynamic shear rheometer (DSR) shown in Figure 3a is the laboratory equipment widely used to characterize the rheological properties of different types of binders and mastics [16]. As shown in Figures 2 and 3b, a disc of the tested material is sandwiched between two plates at the desired temperature. In this test, a sinusoidal force (stress) is applied to the sample, and the deformation (strain) is recorded. The induced strain has a sinusoidal form with a time lag that comes from the material's viscose deformation. This time lag is called the phase angle (δ). The frequency sweep configuration is adopted to obtain the response of fluid-like materials to different loading frequencies by measuring the complex shear modulus and phase angle at different frequencies.  Deformation of Burger's model is the sum of the deformations of Kelvin's elements and Maxwell's elements. This means that the total deformation can be written as: Deformation of the Maxwell element, in turn, comes from the dashpot and the spring: The first and second derivatives of Equation (4) can be written as: . u " . u k`. u mk`. u mc (5) ..
The force-displacement equation of Burger's model can be written as: Although Burger's model is widely used by some researchers, mostly to simulate the mechanical behavior of asphalt mastic, it has not been publicly implemented in YADE [15]. Therefore, authors use an in-house implementation of this model in YADE, creating a cohesive viscoelastic model to simulate pitch. A new type of material, the so-called CohBurgersMat is defined in YADE. This material takes four of Burger's parameters in normal direction and four of Burger's parameters in shear direction. Interactions of elements having CohBurgersMat are governed by Equation (7). Deformation of Burger's model is the sum of the deformations of Kelvin's elements and Maxwell's elements. This means that the total deformation can be written as: Deformation of the Maxwell element, in turn, comes from the dashpot and the spring: The first and second derivatives of Equation (4) can be written as: The force-displacement equation of Burger's model can be written as: Although Burger's model is widely used by some researchers, mostly to simulate the mechanical behavior of asphalt mastic, it has not been publicly implemented in YADE [15]. Therefore, authors use an in-house implementation of this model in YADE, creating a cohesive viscoelastic model to simulate pitch. A new type of material, the so-called CohBurgersMat is defined in YADE. This material takes four of Burger's parameters in normal direction and four of Burger's parameters in shear direction. Interactions of elements having CohBurgersMat are governed by Equation (7).
Viscoelastic properties of a material are normally presented by two parameters of a complex shear modulus, G*, and phase angle, δ [13]. These parameters can be measured using a parallel-plate measuring system as schematically shown in Figure 2. The dynamic shear rheometer (DSR) shown in Figure 3a is the laboratory equipment widely used to characterize the rheological properties of different types of binders and mastics [16]. As shown in Figures 2 and 3b, a disc of the tested material is sandwiched between two plates at the desired temperature. In this test, a sinusoidal force (stress) is applied to the sample, and the deformation (strain) is recorded. The induced strain has a sinusoidal form with a time lag that comes from the material's viscose deformation. This time lag is called the phase angle (δ). The frequency sweep configuration is adopted to obtain the response of fluid-like materials to different loading frequencies by measuring the complex shear modulus and phase angle at different frequencies.  Viscoelastic properties of a material are normally presented by two parameters of a complex shear modulus, G*, and phase angle, δ [13]. These parameters can be measured using a parallel-plate measuring system as schematically shown in Figure 2. The dynamic shear rheometer (DSR) shown in Figure 3a is the laboratory equipment widely used to characterize the rheological properties of different types of binders and mastics [16]. As shown in Figures 2 and 3b, a disc of the tested material is sandwiched between two plates at the desired temperature. In this test, a sinusoidal force (stress) is applied to the sample, and the deformation (strain) is recorded. The induced strain has a sinusoidal form with a time lag that comes from the material's viscose deformation. This time lag is called the phase angle (δ). The frequency sweep configuration is adopted to obtain the response of fluid-like materials to different loading frequencies by measuring the complex shear modulus and phase angle at different frequencies. The following equations are used to calculate G*: in which τ is the maximum shear stress, is the maximum applied torque, is the radius of the specimen, θ is the rotation angle, and ℎ is the thickness of the specimen. Finally, the complex shear modulus is obtained from: * = τ γ The phase angle (δ), in turn, was determined at different frequencies by measuring the delay in seconds between the peaks in the stress and strain functions.
Having obtained * and δ, the storage and loss moduli at 150 °C for each frequency are calculated by:

Experimental
Dynamic shear tests were conducted in strain-controlled mode using an ARES rheometer (see Figure 3). Frequency sweep configuration was adopted at 150 °C from 0.06 Hz to 60 Hz. Complex shear modulus and phase angle were measured at each frequency. Then, the obtained data were used to estimate the four parameters of Burger's model of pitch at 150 °C.

Numerical Method
The open-source discrete element code, YADE [15], was used in this work to simulate the DSR test. The geometry of the numerical model of the DSR test of pitch is shown is Figure 4. Pitch is modeled by an assembly of spheres with the radius of 0.08 mm. Spheres are generated in hexagonal closed pack configuration to make up a disc of 12 mm in diameter and 2 mm in thickness. Size of the elements has been chosen considering the balance between the resolution of the model (size of the elements compared to the size of the sample) and the computation time. The following equations are used to calculate G*: in which τ max is the maximum shear stress, T is the maximum applied torque, r is the radius of the specimen, θ is the rotation angle, and h is the thickness of the specimen. Finally, the complex shear modulus is obtained from: The phase angle (δ), in turn, was determined at different frequencies by measuring the delay in seconds between the peaks in the stress and strain functions.
Having obtained G˚and δ, the storage and loss moduli at 150˝C for each frequency are calculated by:

Experimental
Dynamic shear tests were conducted in strain-controlled mode using an ARES rheometer (see Figure 3). Frequency sweep configuration was adopted at 150˝C from 0.06 Hz to 60 Hz. Complex shear modulus and phase angle were measured at each frequency. Then, the obtained data were used to estimate the four parameters of Burger's model of pitch at 150˝C.

Numerical Method
The open-source discrete element code, YADE [15], was used in this work to simulate the DSR test. The geometry of the numerical model of the DSR test of pitch is shown is Figure 4. Pitch is modeled by an assembly of spheres with the radius of 0.08 mm. Spheres are generated in hexagonal closed pack configuration to make up a disc of 12 mm in diameter and 2 mm in thickness. Size of the elements has been chosen considering the balance between the resolution of the model (size of the elements compared to the size of the sample) and the computation time. A small compressive load was initially applied to ensure the contacts between the plates and the testing material. Frequency sweep dynamic shear tests were run at frequencies ranging from 0.06 Hz to 60 Hz, according to the DSR experiments.
Equations (8)-(12) were then used to calculate the rheological parameters. A cohesive contact model based on Burger's model (Equation (7)) is used to simulate the viscoelastic behavior of pitch at 150 °C. The spherical elements representing pitch have Burger's model in both shear and normal direction in all contacts. In this work, the same values were used for Burger's model of pitch in both shear and normal directions.
In 2011, Liu et al. [17] proposed a new method in using the frequency-temperature superposition (FTS). Traditionally, after conducting dynamic modulus tests at different temperatures, FTS is used in calculating shift factors to build up the master curve at the reference temperature. However, in the approach proposed by Liu et al., * and δ are measured at different frequencies at desired temperatures in the laboratory and the Burger's model parameters are obtained. Then, Burger's dashpot viscosities are modified to predict * and δ at amplified frequencies. This method is of considerable interest for discrete element simulations as for example by using the amplifying coefficient of ξ = 1000, dynamic modulus test at the frequency of 0.06 Hz can be simulated at f = 60 Hz, resulting in significantly reduced computation time. In the present work, ξ = 1000 was used in the DSR tests simulations.

DSR of Pitch and Model Verification
Experimental data of the DSR tests of pitch at 150 °C are given in Figure 5. The obtained data is then used to fit the Burger's model parameters. The fitting procedure is performed by the solver option in Microsoft Excel by minimizing the following function: where ′ and " are respectively the storage and loss moduli experimentally measured at jth frequency; and " are respectively the predicted storage and loss moduli at jth frequency; and m is the number of data points, which is 16.  A small compressive load was initially applied to ensure the contacts between the plates and the testing material. Frequency sweep dynamic shear tests were run at frequencies ranging from 0.06 Hz to 60 Hz, according to the DSR experiments.
Equations (8)-(12) were then used to calculate the rheological parameters. A cohesive contact model based on Burger's model (Equation (7)) is used to simulate the viscoelastic behavior of pitch at 150˝C. The spherical elements representing pitch have Burger's model in both shear and normal direction in all contacts. In this work, the same values were used for Burger's model of pitch in both shear and normal directions.
In 2011, Liu et al. [17] proposed a new method in using the frequency-temperature superposition (FTS). Traditionally, after conducting dynamic modulus tests at different temperatures, FTS is used in calculating shift factors to build up the master curve at the reference temperature. However, in the approach proposed by Liu et al., G˚and δ are measured at different frequencies at desired temperatures in the laboratory and the Burger's model parameters are obtained. Then, Burger's dashpot viscosities are modified to predict G˚and δ at amplified frequencies. This method is of considerable interest for discrete element simulations as for example by using the amplifying coefficient of ξ " 1000, dynamic modulus test at the frequency of 0.06 Hz can be simulated at f = 60 Hz, resulting in significantly reduced computation time. In the present work, ξ " 1000 was used in the DSR tests simulations.

DSR of Pitch and Model Verification
Experimental data of the DSR tests of pitch at 150˝C are given in Figure 5. The obtained data is then used to fit the Burger's model parameters. The fitting procedure is performed by the solver option in Microsoft Excel by minimizing the following function: where G 1 0 j and G 2 0 j are respectively the storage and loss moduli experimentally measured at jth frequency; G 1 j and G 2 j are respectively the predicted storage and loss moduli at jth frequency; and m is the number of data points, which is 16.
where ′ and " are respectively the storage and loss moduli experimentally measured at jth frequency; and " are respectively the predicted storage and loss moduli at jth frequency; and m is the number of data points, which is 16.   Table 1. These parameters are called macroscale parameters, and they determine the global material behavior. The Burger's model parameter values for each contact in the DEM model of pitch, however, depend on the size of two contacting spheres. In our implementation of Burger's model on YADE, these parameters are obtained by the following equation for each contact: in which P M is the macroscale parameter of Burger's model as given in Table 1, L is the diameter of the element, and P m is the microscale parameters of Burger's model of pitch. For a contact formed by two overlapping elements, however, the contact model parameters are obtained as: In the DEM code here, as a contact is formed, its rheological parameters are obtained from Equation (15) and then the force-displacement law, Equation (7), is solved at each cycle as long as the contact is active.
Using the values obtained for pitch and the configuration previously explained in Section 3, DSR tests were simulated in YADE. The DEM model of the DSR test was calibrated by changing the stiffness of the loading plate. Stiffness of the plate was modified to obtain the complex shear modulus of pitch at f = 60 Hz. Then, this value was kept constant in DSR tests at all other frequencies. Using this method, stiffness of the loading plate was determined as 40,000 N/m. Figure 6 presents the simulation results of a DSR test of pitch in terms of stress and strain curves by which, using previously mentioned equations, G˚and δ are calculated. Experimentally measured rheological data (G*, G', and G") of pitch at 150˝C in Figures 7-9 are compared with those obtained by the DEM simulation with the viscoelastic parameters given in Table 1.
of pitch at f = 60 Hz. Then, this value was kept constant in DSR tests at all other frequencies. Using this method, stiffness of the loading plate was determined as 40,000 N/m. Figure 6 presents the simulation results of a DSR test of pitch in terms of stress and strain curves by which, using previously mentioned equations, * and δ are calculated. Experimentally measured rheological data (G * , G', and G") of pitch at 150 °C in Figures 7-9 are compared with those obtained by the DEM simulation with the viscoelastic parameters given in Table 1.   of pitch at f = 60 Hz. Then, this value was kept constant in DSR tests at all other frequencies. Using this method, stiffness of the loading plate was determined as 40,000 N/m. Figure 6 presents the simulation results of a DSR test of pitch in terms of stress and strain curves by which, using previously mentioned equations, * and δ are calculated. Experimentally measured rheological data (G * , G', and G") of pitch at 150 °C in Figures 7-9 are compared with those obtained by the DEM simulation with the viscoelastic parameters given in Table 1.    As can be seen from Figures 7-9, the DEM model of pitch is capable of predicting the complex modulus at a wide range of frequencies with a very good precision. However, storage modulus at high frequencies is over-estimated. This comes from the under-estimation of phase angle at higher frequencies. Similar observation has been reported by Liu and You [17]. They reported that the errors in predictions of the complex modulus and phase angle for the amplified frequencies of less than  As can be seen from Figures 7-9, the DEM model of pitch is capable of predicting the complex modulus at a wide range of frequencies with a very good precision. However, storage modulus at high frequencies is over-estimated. This comes from the under-estimation of phase angle at higher frequencies. Similar observation has been reported by Liu and You [17]. They reported that the errors in predictions of the complex modulus and phase angle for the amplified frequencies of less than As can be seen from Figures 7-9 the DEM model of pitch is capable of predicting the complex modulus at a wide range of frequencies with a very good precision. However, storage modulus at high frequencies is over-estimated. This comes from the under-estimation of phase angle at higher frequencies. Similar observation has been reported by Liu and You [17]. They reported that the errors in predictions of the complex modulus and phase angle for the amplified frequencies of less than 25 kHz is around 4%. However, it increases up to 15% for higher frequencies. Using a smaller time-step and not applying the FTS to avoid its probable inertial effects can improve the predictions of the phase angle.

DEM Simulation of Coke/Pitch Pastes
The verified discrete element model of pitch is then used to predict the viscoelastic properties of pitch/coke mixtures at 150˝C. It should be noted that DSR test measurements on coke/pitch mixtures are not always possible or precise due to the small thickness of the sample disc in the DSR tests and limitations of the setup. Thus, a well-designed DEM model is of considerable practical interest.
Calcined coke has particles with irregular shapes and, as used in anodes, a wide size distribution. For discrete element models in the present work, the method previously developed by the authors [18] was applied to model the coke aggregates by overlapping spheres. Coke particles of the size range of 30 + 50 mesh (0.297-0.595 mm) were modeled as clumps in YADE. An example of a modeled coke particle is shown in Figure 10. Four numerical mixes of coke and pitch with 5, 10, 20, and 30 wt % of coke aggregates were generated. Then, the DSR test was simulated on the samples with the same configuration of the onlypitch sample. The DEM model of coke/pitch mixtures is composed of two phases: spheres in hexagonal closed packed (hcp) configuration representing pitch and clumps representing coke aggregates.
To create the coke/pitch sample models, first, clumps of coke aggregates with the intended size distribution and numbers are randomly placed in a disc with a radius of 6 mm and a height of 2 mm. The model is cycled to remove the possible overlaps between aggregates. At this step, zero gravity is used in the model to prevent the particles from sinking to the bottom of the cylindrical container. Then, using spheres with radii of 0.08 mm, the sample volume is tessellated in an hcp configuration.
Since there are already clumps in the space, there will obviously be some overlaps between the clumps and spheres, which must be managed to have a stable sample.
Using a Python script, big clump-sphere overlaps are removed. The code detects the interactions in which the penetration depth is larger than 0.02 mm and deletes the standalone sphere (pitch element). However, if the overlap value is less than 0.02 mm, the size of the standalone sphere is reduced to 0.04 mm without deleting it. An example of a result of this method is shown in Figure 11. Four numerical mixes of coke and pitch with 5, 10, 20, and 30 wt % of coke aggregates were generated. Then, the DSR test was simulated on the samples with the same configuration of the only-pitch sample. The DEM model of coke/pitch mixtures is composed of two phases: spheres in hexagonal closed packed (hcp) configuration representing pitch and clumps representing coke aggregates.
To create the coke/pitch sample models, first, clumps of coke aggregates with the intended size distribution and numbers are randomly placed in a disc with a radius of 6 mm and a height of 2 mm. The model is cycled to remove the possible overlaps between aggregates. At this step, zero gravity is used in the model to prevent the particles from sinking to the bottom of the cylindrical container. Then, using spheres with radii of 0.08 mm, the sample volume is tessellated in an hcp configuration.
Since there are already clumps in the space, there will obviously be some overlaps between the clumps and spheres, which must be managed to have a stable sample.
Using a Python script, big clump-sphere overlaps are removed. The code detects the interactions in which the penetration depth is larger than 0.02 mm and deletes the standalone sphere (pitch element). However, if the overlap value is less than 0.02 mm, the size of the standalone sphere is reduced to 0.04 mm without deleting it. An example of a result of this method is shown in Figure 11.
Since there are already clumps in the space, there will obviously be some overlaps between the clumps and spheres, which must be managed to have a stable sample.
Using a Python script, big clump-sphere overlaps are removed. The code detects the interactions in which the penetration depth is larger than 0.02 mm and deletes the standalone sphere (pitch element). However, if the overlap value is less than 0.02 mm, the size of the standalone sphere is reduced to 0.04 mm without deleting it. An example of a result of this method is shown in Figure 11.   Figure 11a shows the first step in which only clumps are generated. In Figure 11b, pitch spheres have been added, and the overlapping spheres have been deleted or reduced in diameter to avoid huge overlaps. The resulting paste then undergoes a triaxial compression to have an integrated compacted sample.
In Figure 12, the numerical sample of pitch and coke mixture ready for the DSR tests is shown. The dimensions and test procedure were the same as those for pitch. Four pastes of coke and pitch mixture were numerically prepared with 5, 10, 20, and 30 wt % of coke particles. Size distribution of coke aggregates was the same for all samples in the range of´30 + 50 mesh.
Materials 2016, 9,334 9 of 12 Figure 11a shows the first step in which only clumps are generated. In Figure 11b, pitch spheres have been added, and the overlapping spheres have been deleted or reduced in diameter to avoid huge overlaps. The resulting paste then undergoes a triaxial compression to have an integrated compacted sample.
In Figure 12, the numerical sample of pitch and coke mixture ready for the DSR tests is shown. The dimensions and test procedure were the same as those for pitch. Four pastes of coke and pitch mixture were numerically prepared with 5, 10, 20, and 30 wt % of coke particles. Size distribution of coke aggregates was the same for all samples in the range of −30 + 50 mesh. The effect of the content of coke particles on the complex shear modulus of coke/pitch pastes has been shown in Figure 13. As the content of coke in the mixture increases, the complex shear modulus increases. This raise in * is more pronounced for higher frequencies. However, the addition of coke particles does not have a meaningful effect on the value phase angle and, as the content of the coke particles increases, both storage and loss moduli increase. Similar results have been reported by Pasquino [19] for glass bead suspensions in viscoelastic polymers. The effect of the content of coke particles on the complex shear modulus of coke/pitch pastes has been shown in Figure 13. As the content of coke in the mixture increases, the complex shear modulus increases. This raise in G˚is more pronounced for higher frequencies. However, the addition of coke particles does not have a meaningful effect on the value phase angle and, as the content of the coke particles increases, both storage and loss moduli increase. Similar results have been reported by Pasquino [19] for glass bead suspensions in viscoelastic polymers.
The effect of the content of coke particles on the complex shear modulus of coke/pitch pastes has been shown in Figure 13. As the content of coke in the mixture increases, the complex shear modulus increases. This raise in * is more pronounced for higher frequencies. However, the addition of coke particles does not have a meaningful effect on the value phase angle and, as the content of the coke particles increases, both storage and loss moduli increase. Similar results have been reported by Pasquino [19] for glass bead suspensions in viscoelastic polymers.  Figures 13-15 show that there is a considerable rise in the dynamic moduli of pitch by adding coke particles. In Figure 14, results have also been compared with two equations developed for the case of spherical rigid particles in Newtonian fluids, the Hashin-Shtrikman [20] and Krieger-Dougherty [21] equations, which are respectively written as: where φ is the volume content of particles.
In Figure 14, the ( ) ratio for pitch mixed coke particles are larger than those predicted by the equations. It should be noted that these equations were proposed for Newtonian fluids including  Figures 13-15 show that there is a considerable rise in the dynamic moduli of pitch by adding coke particles. In Figure 14, results have also been compared with two equations developed for the case of spherical rigid particles in Newtonian fluids, the Hashin-Shtrikman [20] and Krieger-Dougherty [21] equations, which are respectively written as: where ϕ is the volume content of particles. In Figure 14, the G1pϕq G1 0 ratio for pitch mixed coke particles are larger than those predicted by the equations. It should be noted that these equations were proposed for Newtonian fluids including rigid spherical particles. However, pitch is a non-Newtonian viscoelastic fluid above its softening point. Added particles here are not spherical either, and they are relatively coarse. Therefore, the case seems complex, and we have found no analytical or empirical equation to describe the effects of adding rigid particles to a viscoelastic fluid. Performing experimental measurement was not possible due to the size of the particles (in the range of 0.297-0.595 mm) and the stiffness of the sample, which exceeded the limits of our setup. Therefore, a mix of pitch with fine particles smaller than 150 micron was made for experimental measurements. The results of DSR tests on these samples compared to the above-mentioned equations are given in Figure 16. Again, it can be seen that the experimental point of the G1 G1 0 ratio for pitch/coke mixtures are larger than the upper limit predicted by the Krieger-Dougherty relation. Thus, it can be concluded that these models may not be appropriate to study the rheological parameters of particles mixed with non-Newtonian viscoelastic fluids. rigid spherical particles. However, pitch is a non-Newtonian viscoelastic fluid above its softening point. Added particles here are not spherical either, and they are relatively coarse. Therefore, the case seems complex, and we have found no analytical or empirical equation to describe the effects of adding rigid particles to a viscoelastic fluid. Performing experimental measurement was not possible due to the size of the particles (in the range of 0.297-0.595 mm) and the stiffness of the sample, which exceeded the limits of our setup. Therefore, a mix of pitch with fine particles smaller than 150 micron was made for experimental measurements. The results of DSR tests on these samples compared to the above-mentioned equations are given in Figure 16. Again, it can be seen that the experimental point of the ratio for pitch/coke mixtures are larger than the upper limit predicted by the Krieger-Dougherty relation. Thus, it can be concluded that these models may not be appropriate to study the rheological parameters of particles mixed with non-Newtonian viscoelastic fluids.

Conclusions
A cohesive viscoelastic contact model based on Burger's model was implemented in YADE. The rheological behavior of pitch was measured at 150˝C using a dynamic shear rheometer. The obtained data were then used to calibrate the DEM model parameters for pitch. DSR test of pitch was then simulated by a three-dimensional DEM model in which pitch is modeled by an assembly of spheres of radius of 0.08 mm in hcp configuration.
Results confirm that Burger's model is a superior choice to describe the complex rheological behavior of pitch. The simulation results are in a very good agreement with the experimental data. However, the storage modulus is over-estimated at high frequencies. It is believed that using a smaller time-step and not applying the frequency-temperature superposition (which of course results in very long computation times) can provide a better prediction of phase angle and thus a better prediction of storage modulus at higher frequencies.
Calibrated DEM model of pitch was then used to predict the rheological properties of coke/pitch pastes with different content of coke aggregates. Since experimental measurements using a DSR machine was not possible, the proposed DEM numerical model is of great practical interest. Results showed that as the content of coke aggregates increases, both storage and loss moduli of the mixture increase, resulting in an near-constant phase angle. The rise in complex modulus value is more pronounced at higher frequencies. Comparison of the obtained result with the available literature shows that the empirical models developed for the case of rigid spheres dispersed in Newtonian fluids fail to accurately predict the effects of adding coke particles to pitch. More investigations are required to elucidate the hydrodynamic interactions between two particles in a viscoelastic suspension and the rheological behavior of the mixture.