# Particle dynamics simulation of wet granulation in a rotating drum

Thanh Trung Vo  Saeid Nezamabadi<EMAIL_ADDRESS>Patrick
Mutabaruka  Jean-Yves Delenne  Edouard Izard  Roland Pellenq  Farhang Radjai
LMGC, Université de Montpellier, CNRS, Montpellier, France.  Bridge and Road
Department, Danang Architecture University, 553000 Da Nang, Vietnam.
$\langle\mbox{MSE}\rangle^{2}$, UMI 3466 CNRS-MIT, MIT Energy Initiative, 77
Massachusetts Avenue, Cambridge 02139, USA.  IATE, UMR1208 INRA - CIRAD -
Université de Montpellier - SupAgro, 34060 Montpellier, France.  ArcelorMittal
R&D Maizières, Voie Romaine, F-57283, Maizières-Lès-Metz, France.

###### Abstract

We simulate the granulation process of solid spherical particles in the
presence of a viscous liquid in a horizontal rotating drum by using molecular
dynamics simulations in three dimensions. The numerical approach accounts for
the cohesive and viscous effects of the binding liquid, which is assumed to be
transported by wet particles and re-distributed homogeneously between wet
particles in contact. We investigate the growth of a single granule introduced
into the granular bed and the cumulative numbers of accreted and eroded
particles as a function of time for a range of values of material parameters
such as mean particle size, size polydispersity, friction coefficient and
liquid viscosity. We find that the granule growth is an exponential function
of time, reflecting the decrease of the number of free wet particles. The
influence of material parameters on the accretion and erosion rates reveals
the nontrivial dynamics of the granulation process. It opens the way to a
granulation model based on realistic determination of particle-scale
mechanisms of granulation.

###### keywords:

granular matter , granulation , capillary bond , Discrete Element Method ,
rotating drum

††journal: Powder Technology

##  1 Introduction

Wet granulation or agglomeration of fine solid particles into larger
particles, called granules or agglomerates, is a widespread technique in
industrial processes such as the manufacture of pharmaceuticals [1],
fertilisers and food products [2, 3], powder metallurgy [4] and iron-making
[5, 6, 7]. The increased size of the granules modifies the rheological
properties of the granular material and may improve flow properties, reduce
the segregation of different types of primary particles or enhance the
permeability for the interstitial gas between grains [3, 8, 9, 10, 11, 12,
13]. The wetting of primary particles is achieved either by mixing them with a
binding liquid prior to the process or by dripping or spraying the liquid to
the material during the process [13, 14, 15, 16, 17, 18, 19, 20, 21]. The
granules nucleate as a result of the collisional-frictional/capillary-viscous
interactions of wet primary particles and increase in size by incorporation of
the available liquid, accretion of primary finer particles and coalescence
with other granules [22]. The existing granules may disappear or keep their
size depending on the amount of available liquid and the rate of erosion as
compared to that of accretion and coalescence [23].

Given the large number of parameters involved in the granulation process, its
detailed physical mechanisms and their relative importance for the resulting
properties of are complex. One may distinguish two different groups of
parameters: process and material [8, 9, 24, 25, 26]. The material parameters
are the properties of the binding liquid and raw material such as liquid
viscosity, primary particle size distribution, mean particle size and friction
coefficient of primary particles [19, 27, 28]. The process parameters are
related to the method of mixing solid particles with the binding liquid and
the corresponding operating parameters. There are different types of
granulators using fluidized beds, high shear or low shear in planetary devices
or rotating drums or disks [8]. Major process parameters in all types of
granulators are the liquid volume [7, 24, 27], granulator size, rotation
speed, inclination angle [29, 30, 31] and filling rate [32]. The granulation
process often needs to be optimized by playing with all these parameters in
order to produce granules of high density, homogeneous distribution of primary
raw particles, a targeted mean size and high strength [33, 34, 35, 36, 37].

The granular dynamics has not the same characteristics in different granulator
geometries and may favor more or less collisional or frictional contacts
between primary particles and influence the redistribution and transport of
the binding liquid. The granulation process is easier to model and control
when the agglomeration is governed by binary collisions between particles, as
in granulators based on fluidized bed or high shearing by impellers. Such
processes have been extensively investigated in application to the
pharmaceutical industry [9]. In contrast, in drum granulators the particles
agglomerate in a downward dense granular flow along inclined rotating drum.
The drum agglomeration has the advantage of being a continuous and robust
process, but since the rheology of dense granular flows is a matter of current
research [38, 39, 40], the agglomeration mechanisms in this geometry remain
quite poorly understood [8]. Granular flows in an inclined rotating drum may
show several flow regimes [41, 42, 43] with the common feature of being dense
and inhomogeneous, and involving inertial effects [14, 18, 23, 44, 45, 46,
47]. A practical difficulty with drum granulation is the in-line monitoring of
the kinetics, making it less amenable to theoretical understanding, which is
required in order to be able to improve drum granulation plants, often
suffering from a significant recycle of undersize and crushed oversize
granules [8].

In this paper, we present extensive simulations of drum granulation by means
of the District Element Method (DEM) [48]. The use of the DEM, in which the
granular material is modeled as an assembly of spherical particles interacting
via frictional/cohesive forces, allows for direct quantification of particle-
scale kinetics and accretion/erosion events. This method has already been
applied to investigate agglomeration in granular shear flows. For example,
Talu et al. [49] introduced a model of wet granulation in which some of the
particles which are assumed to be covered by a binder and therefore sticky
while the rest are dry. The binder-layer between particles dissipates energy
due to viscosity and and allows the particles to stick to another one by the
action of capillary forces. In other reported DEM simulations of the
granulation process, besides capillary and viscous forces, simple empirical
rules are used for progressive wetting of the particles [50, 51, 52]. To
reduce the high computational cost of DEM simulations, some authors have used
the DEM simulation data with a low number of particles to train an artificial
neural network or in conjunction with population balance equations for
application to the large number of particles [51, 53]. The DEM has, however,
never been employed for drum granulation.

In the following, in section 2, we introduce our numerical approach and a
model for the transport and redistribution of the binding liquid. We analyze
in section 3 the evolution of granules as function of the number of drum
rotations by investigating the effects of the process and material parameters
such as rotation speed, Froude number, size ratio between large and small
particles, mean particle size, friction coefficient and viscosity of the
binding liquid, we characterize the dynamics of granulation in terms of the
rate of accretion and erosion. Finally, we conclude in section 4 with a short
summary and possible further research directions.

##  2 Model description and numerical method

Our numerical approach for the simulation of the agglomeration process in a
rotating drum is based on the molecular method and a model for the
redistribution and transport of the binding liquid. We first describe below
the physical assumptions underlying the model. Then, we briefly present the
numerical algorithm with its input parameters and main calculation steps.

###  2.1 Physical assumptions

In the molecular dynamics (MD) method, the equations of motion of all
particles are integrated according to an explicit time-stepping scheme such as
the well-known velocity-Verlet algorithm [48, 54, 55, 56, 57, 58]. A detailed
description of the liquid phase and its interaction with solid particles,
requires sub-particle discretization of the liquid phase and a numerical model
for liquid-gas phase transition [45, 47, 40]. However, such a multi-component
model of partially saturated granular materials is not computationally
efficient for the simulation of the granulation process involving a large
number of primary particles. An efficient alternative approach consists in
accounting for the capillary and lubrication forces between particles as a
well as a particle-scale model for the distribution and transport of the
liquid.

Recent experiments and numerical simulations show that the liquid clusters
condensed from a vapor or introduced by mixing the liquid with grains can be
characterized by their connectivity with the grains. The number of liquid
clusters connected to two grains prevails for low amounts of the liquid. In
this ‘pendular’ state, the liquid is in the form of binary bridges. As the
amount of liquid increases, the clusters involve more and more particles until
a single cluster spans the whole packing. The cohesive effect of the liquid in
thermodynamic equilibrium is controlled by the total wetted surface and the
Laplace pressure. The cohesion rapidly increases as the amount of liquid is
increased in the pendular state, and then it keeps a nearly constant value (or
slightly increases) with increasing amount of the liquid before declining for
large amounts of the liquid [45, 47]. This description assumes, however, that
the particles are in quasi-static equilibrium and the liquid is in
thermodynamic equilibrium. The negative Laplace pressure within the liquid
phase is not uniformly distributed if the system is out of equilibrium.
Furthermore, if the granular material flows (as inside a rotating drum), the
liquid clusters undergo large distortions, and the liquid is continuously re-
distributed as a result of coalescence and separation of liquid clusters. In
practice, a small amount of the added liquid is adsorbed into the particle
rough surfaces and is not directly involved in capillary bonding between
particles.

These features suggest that in a dense granular material, for a broad range of
the amounts of liquid, the cohesive capillary stress is nearly independent of
the amount of liquid, and therefore the effect of liquid volume can be
accounted for by the number of wet particles. In a rotating drum, when the
liquid is poured onto the granular flow, it has not time to diffuse and the
capillary stress leads to the creation of small aggregates of primary
particles that are transported by the granular flow. These “micro-aggregates”
may deform or break up into smaller aggregates. They may also capture more
primary particles or coalesce into larger aggregates if they have an excess
amount of liquid that can be shared with other micro-aggregates as a
consequence of their consolidation under the action of contact forces inside
the granular flow. The initial size of the micro-aggregates is proportional to
that of the droplets but, depending on the wetting method, it can grow rapidly
into granules.

This picture of liquid transport by micro-aggregates, illustrated in Fig. 1,
means that the DEM simulations can be based on the micro-aggregates as wet
units. In its simplest setting, these basic units can be modeled as particles
each transporting a given amount of liquid. Since the number of primary
particles embodied in each micro-aggregate is nearly proportional to the
droplet size, the size distribution of the basic particles can be regarded as
reflecting that of the droplets. When two wet particles (micro-aggregates)
meet, they are subjected to the cohesive action of the spontaneous liquid
bridge appearing between them. In this sense, the interactions between micro-
aggregates are similar to those between particles covered by a liquid layer as
in the method introduced by Talu et al. [49]. In the following, we present in
more detail our numerical method and all model parameters.

Figure 1: Schematic representation of the granulation model.

###  2.2 Numerical method

In the MD method, the particles are modeled as rigid grains interacting via
visco-elastic forces reflecting the contact force and the local contact strain
defined from the relative displacement of particles. Since the particles are
assumed to be rigid, a large repulsive stiffness and hence a high time
resolution are required for the calculation of the interactions between
particles. The motion of each rigid spherical particle $i$ of radius $R_{i}$
is governed by Newton’s second law under the action of normal contact forces
$f_{n}$, tangential contact forces $f_{t}$, capillary forces $f_{c}$, viscous
forces $f_{vis}$ and particle weight $m_{i}\mathbf{g}$ :

$\displaystyle m_{i}\frac{d^{2}\mathbf{r}_{i}}{dt^{2}}$ $\displaystyle=$
$\displaystyle\sum_{j}[(f_{n}^{ij}+f_{c}^{ij}+f_{vis}^{ij})\mathbf{n^{ij}}+f_{t}^{ij}\mathbf{t}^{ij}]+m_{i}\mathbf{g}$
$\displaystyle\mathbf{I}_{i}\frac{d\bm{\omega}_{i}}{dt}$ $\displaystyle=$
$\displaystyle\sum_{j}f_{t}^{ij}\mathbf{c}^{ij}\times\mathbf{t}^{ij}$ (1)

where $\bm{\omega}_{i}$ is the rotation vector of particle $i$, and $m_{i}$,
$\mathbf{I}_{i}$, $\mathbf{r}_{i}$ and $\mathbf{g}$ are the mass, inertia
matrix, position and gravity acceleration vector of particle $i$,
respectively. $\mathbf{n}^{ij}$ denotes the unit vector perpendicular to the
contact plane with particle $j$ and pointing from $j$ to $i$,
$\mathbf{t}^{ij}$ is the unit vector belonging to the contact plane $ij$ and
pointing in the direction opposite to the relative displacement of the two
particles and $\mathbf{c}^{ij}$ is the vector pointing from the center of
particle $i$ to the contact point with particle $j$. The tangential viscous
dissipation, as compared to the normal lubrication force, is neglected [59].
The equations of motion are integrated according to a velocity-Verlet time-
stepping scheme [48, 60].

(a) (b)

Figure 2: (a) Schematic drawing of two different cases of capillary bridges:
particle $i$ in contact with particle $j$ and without contact with particle
$k$; (b) Capillary cohesion force $f_{c}$ as a function of the gap
$\delta_{n}$ between two particles for different values of the liquid volume
$V_{b}$ ($\mu$l) and size ratio $\alpha$.

The normal contact force $f_{n}$ involves two components [61, 62]:

$\displaystyle{f_{n}=f_{n}^{e}+f_{n}^{d}}\ .$ (2)

The normal elastic force $f_{n}^{e}=k_{n}\delta_{n}$ is a linear function of
the normal elastic deflection $\delta_{n}$, where $k_{n}$ is the normal
stiffness constant, and the normal damping force
$f_{n}^{d}=\gamma_{n}\dot{\delta}_{n}$ is proportional to the relative normal
velocity $\dot{\delta}_{n}$, where $\gamma_{n}$ is the normal viscous damping
parameter. These both forces appear only when there is overlap, i.e. for
$\delta_{n}<0$.

As the normal force, the tangential contact force $f_{t}$ is the sum of an
elastic force $f_{t}^{e}=k_{t}\delta_{t}$ and a damping force
$f_{t}^{d}=\gamma_{t}\dot{\delta}_{t}$, where $k_{t}$ is the tangential
stiffness constant, $\gamma_{t}$ denotes the tangential viscous damping
parameter and $\delta_{t}$ and $\dot{\delta}_{t}$ are the tangential
displacement and velocity in contact, respectively. According to the Coulomb
friction law, the contact tangential force is bounded by a force threshold
$\mu f_{n}$, where $\mu$ is the friction coefficient[63, 64, 65, 66]:

$f_{t}=-\\!\min\left\\{(k_{t}\delta_{t}+\gamma_{t}\dot{\delta}_{t}),\mu
f_{n}\right\\}.$ (3)

The capillary attraction force $f_{c}$ between two particles depends on the
liquid volume $V_{b}$ of the liquid bond, liquid-vapor surface tension
$\gamma_{s}$ and particle-liquid-gas contact angle $\theta$; see Fig. 2 [44,
61, 67]. The capillary force is computed from the Laplace-Young equations. In
the pendular state, an approximate solution is given by the expression [63,
68] :

$f_{c}=\left\\{\begin{array}[]{@{}ll@{}}-\kappa\ R,&\mbox{for}\
\delta_{n}<0,\\\ -\kappa\ R\ e^{-\delta_{n}/\lambda},&\mbox{for}\
0\leq\delta_{n}\leq d_{rupt},\\\ 0,&\mbox{for}\
\delta_{n}>d_{rupt},\end{array}\right.$ (4)

where $R=\sqrt{R_{i}R_{j}}$ is the geometrical mean radius of two particles of
radii $R_{i}$ and $R_{j}$ and the capillary force pre-factor $\kappa$ is [63]:

$\kappa=2\pi\gamma_{s}\cos\theta.$ (5)

This force exists up to a debonding distance $d_{rupt}$ given by [61, 62]

$\displaystyle{d_{rupt}=\left(1+\frac{\theta}{2}\right)V_{b}^{1/3}}\ .$ (6)

The characteristic length $\lambda$ is the factor that presents the
exponential falloff of the capillary force in equation (4):

$\lambda=c\ h(r)\ \left(\frac{V_{b}}{R^{\prime}}\right)^{1/2}$ (7)

Here, $R^{\prime}=2R_{i}R_{j}/(R_{i}+R_{j})$ and
$r\\!=\\!\max\\{R_{i}/R_{j};R_{j}/R_{i}\\}$ are the harmonic mean radius and
the size ratio between two particles. The expression (4) nicely fits the
capillary force as obtained from direct integration of the Laplace-Young
equation by setting $h(r)=r^{-1/2}$ and $c\simeq 0.9$ [61, 63, 69]. During the
simulation, the total amount of liquid is evenly re-distributed among all
pairs of wet particles (or micro-aggregates) having a gap below $d_{rupt}$.
The largest value of the capillary force occurs when two particles are in
contact ($\delta_{n}\leq 0$). This is the case of most liquid bonds, and the
capillary force at those contacts is independent of the liquid volume [69].

The normal viscous force $f_{vis}$ is due to the lubrication effect of liquid
bridges between particles. Its classical expression for two smooth spherical
particles is [59, 70]:

$f_{vis}=\frac{3}{2}\pi R^{2}\eta\frac{{v_{n}}}{\delta_{n}},$ (8)

where $\eta$ is the liquid viscosity and $v_{n}$ is the relative normal
velocity assumed to be positive when the gap $\delta_{n}$ is decreasing. This
expression implies that the viscous force diverges when the gap $\delta_{n}$
tends to zero. With this singularity, two rigid particles can not collide in
finite time. However, for slightly rough particles, the surfaces are no longer
parallel and the characteristic size of the asperities allows for collision in
finite time. Hence, we introduce a characteristic length $\delta_{n0}$
representing the size of asperities and assume that the lubrication force is
given by

$f_{vis}=\frac{3}{2}\pi
R^{2}\eta\frac{{v_{n}}}{\delta_{n}+\delta_{n0}}\;\;\;\;\mbox{for}\;\;\delta_{n}>0$
(9)

as long as $\delta_{n}>0$, i.e. for a positive gap. This expression ensures
that the singularity will not occur as long as there is no contact. When
contact occurs, i.e. for $\delta_{n}<0$, we assume that the lubrication force
depends only on the characteristic length, so that

$f_{vis}=\frac{3}{2}\pi
R^{2}\eta\frac{{v_{n}}}{\delta_{n0}}\;\;\;\;\mbox{for}\;\;\delta_{n}\leq 0.$
(10)

In our simulations, we set $\delta_{n0}=5.10^{-4}d_{min}$. This value is small
enough to allow lubrication forces to be effective without leading to its
divergence at contact.

###  2.3 Granulation process model

The rotating drum is a cylinder of length $L$ and diameter $d_{c}$ constructed
geometrically by the juxtaposition of polyhedral rigid elements. Its both ends
are closed by two planes and it can rotate around its axis at given angular
speed $\omega$; see Fig. 3. In all simulations analyzed below, the drum is
horizontal, implying that the gravity is perpendicular to the drum axis.

Figure 3: The geometry of the numerical model of a horizontal rotating drum
used in the granulation simulations.

The drum is filled by allowing 5000 dry particles to fall into it under their
own weights until a stabilized granular bed is obtained. The number of
particles was limited to 5000 in order to simulate a large number of drum
rotations and many runs with different values of the system parameters in a
reasonable computation time. In all simulations, the filling level was
$s_{f}=\frac{2h}{d_{c}}\simeq 0.56$, where $h$ is the filling height of the
granular material inside the drum. We considered three different size classes
in a range $[d_{min},d_{max}]$ with a size ratio $\alpha=d_{max}/d_{min}$.
Each size class has the same total volume so that the size distribution is
uniform in terms of particle volumes. This corresponds to a small number of
large particles and a large number of small particles. This distribution
generally leads to large packing fraction since small particles optimally fill
the pore space between the large particles [71, 72]. The mean particle
diameter $\langle d\rangle$ can also be defined as a function of $\alpha$ by

$\langle d\rangle=d_{min}\frac{2\alpha}{1+\alpha}=d_{max}\frac{2}{1+\alpha}$
(11)

Fig. 4(a) shows an example of the initial state in a system with $\alpha=5$
and $\langle d\rangle$ = 16 $\mu$m. Fig. 4(b) displays the granular flow after
$50$ rotations. At the beginning, all particles are at rest inside the
granulator. With drum rotation, a stable flow configuration is reached after
nearly four rotations.

(a) (b)

Figure 4: Snapshots of the initial state (a) and stable flow (b) in the drum
for $\alpha=5$. The colors show the magnitudes of particle velocities varying
from red (fast particles at the free surface) to blue (slowest particles in
the middle of the drum).

A number $N_{w}$ of particles distributed randomly inside the drum can be
wetted by attributing them a given amount of liquid. These wet particles
collide and nucleate into small granules during the steady flow of the
particles. However, because of the relatively low number of particles, it is
numerically more efficient to start the simulations with a single granule
introduced initially in the center of the granular bed. The granulation
process can then be analyzed by following the particles captured by
(accretion) or subtracted from (erosion) the granule. Obviously, the
coalescence of granules can be investigated only by simulations with a much
larger number of granules. In this paper, we focus on the first option with
granulation around a single wet granule. To define the initial wet granule, we
place a spherical probe in the centre of the granular bed with a radius such
that exactly $N_{g0}=100$ particles are inside the probe. All these particles
are considered to be wet. In addition, we randomly select $N_{w}-N_{g0}=200$
free particles throughout the sample and consider them to be wet. Fig. 5
displays the initial state of the granular bed together with the initial
granule and free wet particles. The capillary and viscous forces are activated
for all wet particles.

Figure 5: Snapshot of the granular bed showing the distribution of dry
particles (in white) and wet particles (in black) both those inside the
initially defined granule in the center of the bed and those randomly
distributed throughout the bed.

The liquid content of wet particles $w=V_{\ell}/V_{g}$, where $V_{\ell}$ is
the amount of the liquid and $V_{g}$ is the particle volume, is assumed to
have the same value for all wet particles. In all our granulation simulations,
we set $w=0.09$, which is sufficient to create wet granules in a horizontal
rotating drum [6]. When two wet particles meet, the volume of the liquid
bridge is calculated from $V_{\ell}=wV_{g}$ transported by each of the
particles. We also assume that there is no excess liquid so that a wet
particle can not form a liquid bridge with a dry particle.

We performed a large number of simulations with different values of $d_{min}$
and $d_{max}$ in the range $[10,1500]$ $\mu$m, different values of friction
coefficient $\mu$ in the range $[0.1,0.9]$ and four values of liquid viscosity
$\eta$. For each set of values, several independent granular beds were
generated and subjected to the granulation process. All values of the system
parameters used in our simulations are listed in Table 1. The choice of most
parameter values is guided either by numerical efficiency or by reference to
the granulation process of iron ores in a rotating drum, e.g. for the density
of particles, contact angle and filling level. However, for data analysis we
rely on dimensionless parameters. For rotating drum a relevant dimensionless
parameter is the Froude number $Fr$ [42, 43, 73]:

$Fr=\frac{\omega^{2}d_{c}}{2g}.$ (12)

The flow regime depends on the value of the Froude number. For dry particles,
the value $Fr=0.5$ leads to a flow regime intermediate between rolling and
cascading regimes [25, 59]. Since these regimes have been established for dry
particles, we investigated the effect of the wet particles by comparing the
velocity profiles for dry particles, on the one hand, and in the presence of
wet particles, on the other hand, for $Fr=0.5$. We observe practically the
same velocity profiles in both cases. Fig. 6 shows the mean particle velocity
$v$ and free surface inclination $\beta$ for the two simulations with
$\alpha=5$. We see that both $v$ and $\beta$ have the same value in both
flows. The only difference is that there are more fluctuations in the presence
of free wet particles that can be attributed to the higher inhomogeneity of
the flow in this case.

Figure 6: The mean velocity $v$ of drum flow for dry particles and in the presence of wet particles with $Fr=0.5$ as a function of the number of rotations of the drum. The inset shows the evolution of the free surface inclination angle $\beta$ as a function of the number of rotations. Table 1: Simulation parameters Parameter | Symbol | Value and Unit
---|---|---
Particle diameters | $d$ | [10,1500] $\mu$m
Density of particles | $\rho$ | 3500 kg.m-3
Size ratios | $\alpha$ | [1,5]
Number of particles | $N_{p}$ | 5000
Filling level | $s_{f}$ | 0.56
Friction coefficient | $\mu$ | [0.1,0.9]
Normal stiffness | $k_{n}$ | 100 N/m
Tangential stiffness | $k_{t}$ | 80 N/m
Normal damping | $\gamma_{n}$ |  $5.10^{-5}$ Ns/m
Tangential damping | $\gamma_{t}$ |  $5.10^{-5}$ Ns/m
Surface tension | $\gamma_{s}$ | 0.021 N/m
Contact angle | $\theta$ | 4.0 degree
Liquid viscosity | $\eta$ | [10;20;40;60] mPa.s
Time step | $\delta t$ | 10-7 s

##  3 Parametric study of granulation

In this section, we are interested in the evolution of the granule size, in
terms of the total number of particles $N_{g}$ embodied in the granule, as
well as the relative contributions of the accretion and erosion events. The
cumulative accretion is the number $N_{g}^{+}$ of free wet particles captured
by the granule whereas the cumulative erosion $N_{g}^{-}$ is the number of wet
particles leaving the granule. The rates of these events depend on the
relative importance of force chains and cohesive stresses acting on the
granule. The values of the process parameters affect the rates so that the
granule may grow at different rates. When the rate becomes negative, the
granule initially inserted into the granular bed will disappear by excess
erosion. Unless explicitly stated, the liquid properties are those of water
($\eta=1$ mPa.s and $\gamma_{s}=0.072$ N/m). From the parametric study, we
will determine the phase diagram of granule growth for polydispersity vs. mean
particle size.

###  3.1 Growth, accretion and erosion

Figure 7 shows the evolution of the granule size as a function of the number
of drum rotations $N_{r}=\omega t/2\pi$, where $t$ is the granulation time,
for different values of the size ratio $\alpha$, $Fr=0.5$ and $\mu=0.5$. In
these simulations, $d_{min}$ is kept constant and equal to 10 $\mu$m. This
means that $\alpha$ is increased here by increasing $d_{max}=\alpha d_{min}$.
The granule grows almost exponentially with $N_{r}$ at a rate that increases
with $\alpha$. In other words, the increased size polydispersity enhances
granulation. This effect is more spectacular when compared to the mono-
disperse case ($\alpha=1$) in which the granule growth is negligibly small
after 50 rotations. As we are interested in this paper in the evolution of
captured and eroded particles, the granule size $N_{g}$ is expressed in terms
of the number of wet particles in the granule. Even for broad size
polydispersity, we find that the total volume of wet particles belonging to
the granule is a linear function of their number, as shown in Fig. 8, up to a
factor that depends on the size polydispersity. Hence, the trends of $N_{g}$
investigated below are similar if the granule size is measured in terms of the
total volume of the granule or the total volume of particles in the granule.

The cumulative accreted and eroded particles are plotted in Fig. 9 as a
function of $N_{r}$ only for poly-disperse samples in which we have a
significant number of particles captured and eroded. The accretion grows
exponentially at a negative rate whereas erosion is a linear function. We will
see below with other values of material parameters that the erosion can grow
in a non-linear way and faster than accretion so that the linear evolution
observed in Fig. 9 may be considered as a first-order effect in the limit of
low erosion rates. Fig. 9 also shows that both accretion and erosion increase
with $\alpha$ to an extent that is higher for accretion than for erosion.

Figure 7: Evolution of the granule size $N_{g}$ (in number of particles) for
different values of size ratio $\alpha$ and $d_{min}=10$ $\mu$m. The solid
lines are exponential fits given by equation (16). Figure 8: Variation $\Delta
V$ of the volume of particles inside the granule as a function of their number
$N_{g}$ for different values of size ratio $\alpha$ and $d_{min}=10$ $\mu$m.

(a) (b)

Figure 9: Cumulative accretion (a) and cumulative erosion (b) of particles for
different values of size ratio $\alpha$. The lines are fitting forms given by
equations (17) and (18).

The exponential increase of the number of accreted particles $N_{g}^{+}$ is a
consequence of the decreasing number of available free wet particles in the
granular bed while they are captured by the granule. The number of free wet
particles is given by $N_{w}-N_{g}$, where $N_{w}$ is the total number of wet
particles including those belonging to the granule. This is equivalent to the
decrease of the available liquid for granulation as the granules grow. Hence,
in the steady flow, we may assume that the variation $\Delta N_{g}^{+}$ of the
captured particles is proportional to the current number $N_{w}-N_{g}$ of wet
particles and to the elapsed time $\Delta t$ or angular rotation $\Delta
N_{r}=\omega\Delta t$:

$\Delta N_{g}^{+}=k^{+}\frac{N_{w}-N_{g}}{N_{w}}\Delta N_{r},$ (13)

where $k^{+}$ is the relative accretion rate. As to the number $N_{g}^{-}$ of
eroded particles, we assume that its rate $k^{-}$ is constant :

$\Delta N_{g}^{-}=k^{-}\Delta N_{r},$ (14)

These equations are consistent with our numerical data shown in Fig. 9
although we expected the number of eroded particles to be proportional to the
number of particles at the surface of the granule. This effect may reflect the
fact that the average curvature of the granule surface declines as its size
increases so that the particles lying at the surface of the granule are more
strongly attached to the granule and/or less subjected to the eroding action
of granular flow. This effect may counter-balance the increase of the granule
surface area. However, for much larger granules in number of primary particles
this effect may disappear.

With the above assumptions, the rate equation for the granule size is simply
$\Delta N_{g}=\Delta N_{g}^{+}-\Delta
N_{g}^{-}=(-k^{+}N_{g}/N_{w}+k^{+}-k^{-})\Delta N_{r}$, which leads to a
simple differential equation:

$\frac{dN_{g}}{dN_{r}}=k^{+}\left(1-N_{g}/N_{w}\right)-k^{-}$ (15)

with the following solution:

$N_{g}(N_{r})=N_{g0}+\\{N_{w}(1-k^{-}/k^{+})-N_{g0}\\}(1-e^{\frac{-k^{+}}{N_{w}}N_{r}})$
(16)

This model predicts an exponential growth and an asymptotic granule size
$N_{g}(t\to\infty)=N_{w}(1-k^{-}/k^{+})$, as observed in our simulations. The
steady granulation state corresponds to the condition $\Delta N_{g}^{+}=\Delta
N_{g}^{-}$. For $k^{-}\ll k^{+}$, the final granule embodies nearly all wet
particles: $N_{g}\simeq N_{w}$. On the other hand, the granule disappears if
$k^{-}>k^{+}$.

Figure 10: Fitted values of erosion rate $k^{-}$ and accretion rate $k^{+}$ as
a function of polydispersity parameter $\alpha$.

The evolution of $N_{g}^{+}$ and $N_{g}^{-}$ as a function of drum rotation
can be obtained from equations (16), (13) and (14). We get

$N_{g}^{-}=k^{-}N_{r}$ (17)

and

$N_{g}^{+}=\\{N_{w}(1-k^{-}/k^{+})-N_{g0}\\}(1-e^{\frac{-k^{+}}{N_{w}}N_{r}})+k^{-}N_{r}$
(18)

In view of the present model, the influence of particle size ratio $\alpha$
observed in Figs. 7 and 9 can be interpreted in terms of the accretion and
erosion rates. Fig.10 shows the fitted values of $k^{-}$ and $k^{+}$ as a
function of $\alpha$. The increase of accretion rate $k^{+}$ with $\alpha$ is
rather counter-intuitive since the cohesive strength is inversely proportional
to the mean particle size, which increases here with $\alpha$. This means that
the higher polydispersity, allowing for a better filling of the pore space and
thus higher density of the granule, over-compensates the decrease of the
cohesive stress. But the latter explains the increase of the erosion rate
$k^{-}$, which is quite small compared to $k^{+}$.

###  3.2 Effects of material parameters

We now consider the effect of the mean particle size $\langle d\rangle$, which
directly controls the cohesive stress of wet particles. Fig. 11 displays the
cumulative number of particles for both accretion and erosion in the case
$\alpha=5$, $\mu=0.5$. Note that, in these simulations, the higher values of
$\langle d\rangle$ imply higher values of both $d_{min}$ and $d_{max}$. But,
according to equation (11), $d_{max}$ increases faster than $d_{min}$. Fig. 11
shows that accretion $N_{g}^{+}$ increases as an exponential function of the
number of drum rotations whereas erosion $N_{g}^{-}$ is quasi-linear. As
expected, since the cohesive stress declines, accretion decreases and erosion
increases with increasing $\langle d\rangle$. For $\langle d\rangle=1666$
$\mu$m, erosion is high enough to cancel the effect of accretion, and thus the
granule disappears after 8 rotations. Here, the cumulative erosion does not
grow linearly with time and for this reason the erosion rate is not constant.
For the linear part of erosion, we have $k^{-}\simeq 25$ whereas $k^{+}\simeq
6$. Clearly, $k^{-}$ is strongly dependent on the cohesive stress of the
granule, which declines in inverse proportion to the wet particle mean size.
The cohesive stress of wet particles affects, albeit to a lesser extent, the
accretion rate $k^{+}$.

Figure 11: Evolution of the cumulative number of wet particles for accretion
(solid lines) and erosion (dashed lines) for four different values of the mean
particle diameter $\langle d\rangle$ and size ratio $\alpha=5$, as a function
of the number of drum rotations.

The friction coefficient $\mu$ between particles is a major parameter for
granular flows. Its joint effect with cohesive forces can thus influence the
granulation process. Fig. 12 shows cumulative accretion and erosion for an
increasing value of $\mu$ and exponential fits. Fig. 13 displays $k^{-}$ and
$k^{+}$ as a function of $\mu$. We observe a slight increase of both rates.
The values of $N_{g}^{+}$ and $N_{g}^{-}$ after 50 rotations show that
accretion increases slightly with $\mu$, which may be understood as a
consequence of enhanced capturing of free wet particles when they touch the
granule. Surprisingly, however, $N_{g}^{-}$ increases to a larger extent with
$\mu$ so that the granule growth is slower at larger values of $\mu$. As
erosion is a consequence of interactions between dry particles and boundary
wet particles of the granule, the increase of erosion with $\mu$ may be
understood as an increase of shear forces acting on the boundary particles by
shear flow of dry particles. This suggests that the granulation of rough
particles (with higher friction coefficient) is less efficient than rounded
particles and it should consume more energy.

Figure 12: Cumulative accretion (solid lines) and erosion (dashed lines) of
wet particles for five different values of the friction coefficient $\mu$,
size ratio $\alpha=5$ and $d_{min}=10$ $\mu$m, as a function of the number of
drum rotations. The lines are fitting forms given by equations (17) and (18).
Figure 13: Fitted values of erosion rate $k^{-}$ and accretion rate $k^{+}$ as
a function of friction coefficient $\mu$ for size ratio $\alpha=5$.

This brings us to the effect of liquid viscosity expressed as a lubrication
force between wet particles. We performed several granulation simulations for
different values of $\eta$ in the range $[10,60]$ mPa.s and with same value of
the liquid-vapor surface tension $\gamma_{s}$=21 mN/m. Fig. 14 shows the
evolution of cumulative accretion and erosion for $\alpha=5$. Again, we
observe the exponential increase of $N_{g}^{+}$ vs. the nearly linear increase
of $N_{g}^{-}$ with the number of drum rotations. The accretion increases
slightly with $\eta$ whereas the erosion declines. The decrease of $N_{g}^{-}$
is exactly the opposite effect of friction coefficient in Fig. 12 which causes
an increase of $N_{g}^{-}$. This means that lubrication forces tend to reduce
the shearing effect of the flow on the granule, leading to smaller erosion. On
the other hand, the increase of accretion can be attributed to viscous
dissipation that can enhance the capture of free wet particles by the granule.
Fig. 15 shows $k^{-}$ and $k^{+}$ as a function of $\eta$. The variations of
the rates with $\eta$ are small but their opposite effects tend to enhance
granule growth.

Figure 14: Cumulative accretion (solid lines) and erosion (dashed lines) for
four different values of the liquid viscosity $\eta$ for size ratio $\alpha=5$
and $d_{min}=10$ $\mu$m as a function of the number of drum rotations. The
solid and dashed lines are exponential and linear fits to the data points,
respectively. Figure 15: Erosion rate $k^{-}$ and accretion rate $k^{+}$ as a
function of liquid viscosity $\eta$ for size ratio $\alpha=5$ and $d_{min}=10$
$\mu$m.

###  3.3 Phase diagram

Figure 16: Phase diagram of granule growth in the parametric space of $\langle
d\rangle$ vs. $\alpha$ for $\eta$ = 1 mPa.s and $\mu=0.5$. The granule grows
only in the light blue region and disappears otherwise.

In the simulations reported in this paper, we have a single granule of size
$N_{g0}$ that can increase or decrease in size depending on the values of
various material parameters. In the last section, we analyzed the influence of
several parameters on the accretion and erosion rates. The most crucial issue,
however, in this single-granule problem is the ranges of the values of those
parameters for which the granule will survive and grow, i.e. the phase diagram
of granulation. We consider here only the phase diagram in the parameter space
of the mean particle size $\langle d\rangle$ vs. size ratio $\alpha$. In these
simulations, we set $\gamma_{s}$=72 mN/m and $\mu=0.5$.

Figure 16 displays all the simulated points as a grid in the parameter space
$[\langle d\rangle,\alpha]$ and an approximate frontier between the range of
values for which the granule survives, i.e. granulation is possible, and the
values for which no granule can grow. The liquid viscosity is set to $\eta=1$
mPa.s. We see that for $\langle d\rangle>1200$ $\mu$m no granulation occurs
even for larger values of $\alpha$ (larger polydispersity). For smaller values
of $\langle d\rangle$, the range of the values of $\alpha$ for granulation
increases. Note that the limit on the value of $\alpha$ is related to the fact
that for given $\langle d\rangle$, the increase of $\alpha$ requires the
increase of $d_{max}$ and decrease of $d_{min}$. The increase of $d_{max}$
compensates to some extent the effect of decreasing $d_{min}$ on the cohesive
stress.

##  4 Conclusions

In this paper, we used a 3D molecular dynamics algorithm with a capillary
cohesion law enhanced by the viscous effect of the binding liquid in order to
investigate the granulation process of granular materials in a horizontal drum
rotating about its axis. The system was numerically prepared by pouring solid
particles into a rigid drum composed of polyhedral elements. We considered in
detail the evolution of a single granule introduced from the beginning into
the granular bed for a broad range of parameter values. The liquid was assumed
to be transported by the particles and its amount was defined by the number of
wet free particles randomly distributed inside the granular bed.

We showed that the granule grows almost exponentially with time (or the number
of drum rotations) as a result of a gradual capture of free wet particles by
the granule. The accretion of free wet particles is nearly always an
exponential function of time whereas the number of eroded particles grows
linearly with time. A simple model based on constant accretion and erosion
rates was introduced, predicting the observed exponential increase of granule
size leveling off to a constant value at long times. We investigated the
effects of size ratio, mean particle size, friction coefficient and liquid
viscosity on accretion and erosion of particles. Both accretion and erosion
increase when the size ratio or friction coefficient are increased. Accretion
declines whereas erosion increases when the mean particle size is increased.
Accretion increases whereas erosion declines when liquid viscosity is
increased. We determined the phase diagram of granulation by varying
systematically the size ratio and mean particle size. This will be extended to
other parameters in the future.

Our results are based on a simple system with a relatively small number of
spherical particles. We considered the evolution of a single granule and this
allowed us to perform long-time simulations for a range of values of material
parameters. Clearly, these simulations can be extended to obtain the combined
effects of parameters in phase-space diagrams. The effects of process
parameters such as filling rate, rotation speed and wetting procedure can be
studied, too. Further data analysis is also necessary in order to investigate
the granule consolidation and the problem of nucleation from free wet
particles when no granule is initially introduced. Despite its higher
computational cost, a larger number of particles may allow for the simulation
of multi-granule systems and coalescence phenomena.

## Acknowledgments

We gratefully acknowledge financial support by the Ministry of Education and
Training in Vietnam and Campus France. This work was initiated in
collaboration with ArcelorMittal.

## References

  * [1] S. H. Chien, G. Carmona, L. I. Prochnow, E. R. Austin, Cadmium availability from granulated and bulk-blended phosphate-potassium fertilizers, Journal of Environmental Quality 32 (2003) 1911–1914.
  * [2] K. Saleh, P. Guigon, Coating and encapsulation processes in powder technology, Handbook of Powder Technology, Granulation 11 (2007) 323–375.
  * [3] E. Rondet, M. Delalonde, T. Ruiz, J. P. Desfoursb, Fractal formation description of agglomeration in low shear mixer, Chemical Engineering Journal 164 (1971) 212–219.
  * [4] A. Nosrati, J. Addai-Mensah, D. J. Robinson, Drum agglomeration behavior of nickel laterite ore: Effect of process variables, Hydrometallurgy 125–126 (2012) 90–99.
  * [5] K. V. Sastry, P. Dontula, C. Hosten, Investigation of the layering mechanism of agglomerate growth during drum pelletization, Powder Technology 130 (2003) 231–237.
  * [6] R. Aguado, S. Roudier, L. Delagado, Best available techniques (BAT) reference document for iron and steel production, Joint Research Centre of the European Commission, Luxembourg: Publications Office of the European Union (2013).
  * [7] G. M. Walker, Chapter 4 Drum Granulation Processes, Handbook of Powder Technology 11 (2007) 219 – 254, granulation.
  * [8] K. Saleh, P. Guigon, Handbook of Powder Technology, Granulation, Vol. 11, Elsevier, Amsterdam, 2007, Ch. Coating and encapsulation processes in powder technology, pp. 323–375.
  * [9] P. Suresh, I. Sreedhar, R. Vaidhiswaran, A. Venugopal, Acomprehensive review on process and engineering aspects of pharmaceutical wet granulation, Chemical Engineering Journal 328 (2017) 785–815.
  * [10] S. M. Iveson, J. D. Litster, K. Hapgood, B. J. Ennis, Nucleation, growth and breakage phenomena in agitated wet granulation processes: a review, Powder Technology 117 (1) (2001) 3–39.
  * [11] S. Herminghaus, Dynamics of wet granular matter, Advances in Physics 54 (3) (2005) 221–261.
  * [12] J. Litster, B. Ennis, The Science and Engineering of Granulation Processes, Springer Netherlands 15 (2014).
  * [13] C. Liao, S. Hsiau, S. Wen, Effect of adding a small amount of liquid on density-induced wet granular segregation in a rotating drum, Advanced Powder Technology 27 (4) (2016) 1265–1271.
  * [14] S. Iveson, J. Litster, Growth regime map for liquid-bound granules, AIChE journal 44 (7) (1998) 1510–1518.
  * [15] X. Xiao, Y. Tan, H. Zhang, R. Deng, S. Jiang, Experimental and dem studies on the particle mixing performance in rotating drums: Effect of area ratio, Powder Technology 314 (2017) 182–194.
  * [16] B. Alchikh-Sulaiman, M. Alian, F. Ein-Mozaffari, A. Lohi, S. R.Upreti, Using the discrete element method to assess the mixing of polydisperse solid particles in a rotary drum, Particuology 25 (2016) 133–142.
  * [17] A. Bouwman, M. Henstra, D. Westerman, J. Chung, Z. Zhang, A. In- gram, J. Seville, H. Frijlink, The effect of the amount of binder liquid on the granulation mechanisms and structure of microcrystalline cellulose granules prepared by high shear granulation, International Journal of Pharmaceutics 290 (1) (2005) 129–136.
  * [18] J. Degre‘ve, J. Baeyens, M. V. de Velden, S. D. Laet, Spray-agglomeration of npk-fertilizer in a rotating drum granulator, Powder Technology 163 (3) (2006) 188–195.
  * [19] J. D. Osborne, R. P. Sochon, J. J. Cartwright, D. G. Doughty, M. J. Hounslow, A. D. Salman, Binder addition methods and binder distribution in high shear and fluidised bed granulation, Chemical Engineering Research and Design 89 (5) (2011) 553–559.
  * [20] S. K. Pawar, F. Henrikson, G. Finotello, J. T. Padding, N. G. Deen, A. Jongsma, F. Innings, J. H. Kuipers, An experimental study of dropletparticle collisions, Powder Technology 300 (2016) 157–163.
  * [21] R. Pashminehazar, A. Kharaghani, E. Tsotsas, Three dimensional characterization of morphology and internal structure of soft material agglomerates produced in spray fluidized bed by x-ray tomography, Powder Technology 300 (2016) 46–60.
  * [22] B. J. Ennis, G. Tardos, R. Pfeffer, A microlevel-based characterization of granulation phenomena, Powder Technology 65 (1) (1991) 257–272.
  * [23] F. Stepanek, P. Rajniak, C. Mancinelli, R. Chern, R. Ramachandran, Distribution and accessibility of binder in wet granules, Powder Technology 189 (2) (2009) 376–384.
  * [24] N. Rahmanian, A. Naji, M. Ghadiri, Effects of process parameters on granules properties produced in a high shear granulator, Chemical Engineering Research and Design 89 (5) (2011) 512–518.
  * [25] B. J. Ennis, Agglomeration and size enlargement session summary paper, Powder Technology 88 (3) (1996) 203–225.
  * [26] M. Butensky, D. Hyman, Rotary drum granulation. an experimental study of the factors affecting granule size, Ind. Eng. Chem. Fundam. 10 (1971) 212–219.
  * [27] K. E. Ileleji, Y. Li, R. K. Ambrose, P. H. Doane, Experimental investigations towards understanding important parameters in wet drum granulation of corn stover biomass, Powder Technology 300 (2016) 126–135.
  * [28] R. Ramachandran, J. M. H. Poon, C. F. Sanders, T. Glaser, C. D. Immanuel, F. J. Doyle, J. D. Litster, F. Stepanek, F.-Y. Wang, I. T. Cameron, Experimental studies on distributions of granule size, binder content and porosity in batch drum granulation: Inferences on process modeling requirements and process sensitivities, Powder Technology 188 (2) (2008) 89–101.
  * [29] J. P. Pan, T. J. Wang, J. J. Yao, Y. Jin, Granule transport and mean residence time in horizontal drum with inclined flights, Powder Technology 162 (1) (2006) 50–58.
  * [30] R. Spurling, J. Davidson, D. Scott, The transient response of granular flows in an inclined rotating cylinder, Chemical Engineering Research and Design 79 (1) (2001) 51–61.
  * [31] F. Wang, I. Cameron, A multi-form modeling approach to the dynamics and control of drum granulation processes, Powder Technology 179 (1) (2007) 2–11.
  * [32] J. M. N. T. Gray, Granular flow in partially filled slowly rotating drums, Journal of Fluid Mechanics 441 (2001) 1–29.
  * [33] C. Thornton, Z. Ning, A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres, Powder Technology 99 (2) (1998) 154–162.
  * [34] L. X. Liu, J. D. Litster, S. M. Iveson, B. J. Ennis, Coalescence of de- formable granules in wet granulation processes, AIChE Journal 46 (2000) 529–539.
  * [35] S. Iveson, J. Beathe, N. Page, The dynamic strength of partially saturated powder compacts: the effect of liquid properties, Powder Technology 127 (2002) 149–161.
  * [36] M. Ghadiri, A. D. Salman, M. Hounslow, A. Hassanpour, D. W. York, Editorial: Special issue ? agglomeration, Chemical Engineering Research and Design 89 (5) (2011) 499.
  * [37] J. Rojek, S. Nosewicz, M. Mazdziarz, P. Kowalczyk, K. Wawrzyk, D. Lumelskyj, Modeling of a sintering process at various scales, Procedia Engineering 177 (2017) 263–270.
  * [38] GDR-MiDi, On dense granular flows, Eur. Phys. J. E 14 (2004) 341–365.
  * [39] K. Kamrin, D. L. Henann, Nonlocal modeling of granular flows down inclines, Soft Matter, 2015 11 (2015) 179–185.
  * [40] L. Amarsid, J.-Y. Delenne, P. Mutabaruka, Y. Monerie, F. Perales, F. Radjai, Viscoinertial regime of immersed granular flows, Phys. Rev. E 96 (2017) 012901.
  * [41] I. Govender, Granular flows in rotating drums: A rheological perspective, Minerals Engineering 92 (2016) 168–175.
  * [42] R. Yang, A. Yu, L. McElroy, J. Bao, Numerical simulation of particle dynamics in different flow regimes in a rotating drum, Powder Technology 188 (2) (2008) 170–177.
  * [43] J. Mellmann, The transverse motion of solids in rotating cylinders?forms of motion and transition behavior, Powder technology 118 (3) (2001) 251–270.
  * [44] G. Lian, C. Thornton, M. Adams, A theoretical study of the liquid bridge forces between two rigid spherical bodies, Journal of Colloid and Interface Science 161 (1993) 138–147.
  * [45] M. Scheel, R. Seemann, M. Brinkmann, M. D. Michiel, A. Sheppard, S. Herminghaus, Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime, Journal of Physics: Condensed Matter 20 (49) (2008) 494236.
  * [46] C. Willett, M. Adans, S. Johnson, J. Seville, Capillary bridges between two spherical bodies, Langmuir 16 (2000) 9396–9405.
  * [47] J.-Y. Delenne, V. Richefeu, F. Radjai, Liquid clustering and capillary pressure in granular media, Journal of Fluid Mechanics 762.
  * [48] F. Radjai, F. Dubois, Discrete-element modeling of granular materials, Wiley-Iste, 2011.
  * [49] I. Talu, G. I. Tardos, M. I. Khan, Computer simulation of wet granulation, Powder Technology 110 (2000) 59–75.
  * [50] E. L. Chan, K. Washino, G. K. Reynolds, B. Gururajan, M. J. Hounslow, A. D. Salman, Blade-granule bed stress in a cylindrical high shear granulator: Further characterisation using dem, Powder Technology 300 (2016) 92–106.
  * [51] P. Lau, M. Kind, Cfd-pbe simulation to predict particle growth in a fluidized bed melt granulation batch process, Powder Technology 300 (2016) 28–36.
  * [52] S. Sarkar, B. Chaudhuri, Dem modeling of high shear wet granulation of a simple system, Asian Journal of Pharmaceutical Sciences 13 (2018) 220–228.
  * [53] D. Barrasso, A. Tamrakar, R. Ramachandran, A reduced order pbm?ann model of a multi-scale pbm?dem description of a wet granulation process, Chemical Engineering Science 119 (2014) 319–329.
  * [54] V. Richefeu, M. S. El Youssoufi, E. Aze?ma, F. Radjai, Force transmission in dry and wet granular media, Powder Technology 190 (2009) 258–263.
  * [55] H. Matuttis, S. Luding, H. Herrmann, Discrete element simulations of dense packings and heaps made of spherical and non-spherical particles, Powder Technology 109 (1) (2000) 278–292.
  * [56] C. Thornton, Numerical simulations of deviatoric shear deformation of granular media, Geotechnique 50 (1) (2000) 43–53.
  * [57] F. Radjai, I. Preechawuttipong, R. Peyroux, Cohesive granular texture, in: P. Vermeer, S. Diebels, W. Ehlers, H. Herrmann, S. Luding, E. Ramm (Eds.), Continuous and discontinuous modeling of cohesive frictional materials, Springer Verlag, Berlin, 2001, pp. 148–159.
  * [58] H. J. Herrmann, S. Luding, Modeling granular media with the computer, Continuum Mechanics and Thermodynamics 10 (1998) 189–231.
  * [59] G. Lefebvre, P. Jop, Erosion dynamics of a wet granular medium, Physical Review E : Statistical, Nonlinear, and Soft Matter Physics 8 (2013) 032205.
  * [60] J. Duran, A. Reisinger, P. de Gennes, Sands, Powders, and Grains: An Introduction to the Physics of Granular Materials, Partially Ordered Systems, Springer New York, 1999.
  * [61] V. Richefeu, M. El Youssoufi, F. Radjai, F., Shear strength properties of wet granular materials, Physical Review E 73 (2006) 051304.
  * [62] T-Trung. Vo, P. Mutabaruka, J-Y. Delenne, S. Nezamabadi, F. Radjai, Strength of wet agglomerates of spherical particles: effects of friction and size distribution, EPJ Web Conf. 140 (2017) 08021.
  * [63] V. Richefeu, F. Radjai, M. S. E. Youssoufi, Stress transmission in wet granular materials., Eur. Phys. J. E 21 (2007) 359–369.
  * [64] J. Scha?fer, S. Dippel, D. E. Wolf, Force schemes in simulations of granular materials, J. Phys. I France 6 (1996) 5–20.
  * [65] S. Dippel, G. G. Batrouni, D. E. Wolf, How transversal fluctuations affect the friction of a particle on a rough incline, Phys. Rev. E 56 (1997) 3645–3656.
  * [66] S. Luding, Collisions and contacts between two particles, in: H. J. Herrmann, J.-P. Hovi, S. Luding (Eds.), Physics of dry granular media - NATO ASI Series E350, Kluwer Academic Publishers, Dordrecht, 1998, p. 285.
  * [67] V. Richefeu, M. S. El Youssoufi, F. Radja??, Shear strength of unsaturated soils: Experiments, dem simulations, and micromechanical analysis, in: Theoretical and Numerical Unsaturated Soil Mechanics, Springer, 2007, pp. 83–91.
  * [68] T-Trung. Vo, P. Mutabaruka, S. Nezamabadi, J.-Y. Delenne, E. Izard, R. Pellenq, F. Radjai, Mechanical strength of wet particle agglomerates, Mechanics Research Communications 92 (2018) 1–7.
  * [69] F. Radjai, V. Richefeu, Bond anisotropy and cohesion of wet granular materials, Philosophical Transactions of the Royal Society A 367 (2009) 5123–5138.
  * [70] J. Happel, H. Brenner, Low Reynolds Number Hydrodynamics, Martinus Nijhoff Publishers, 1983.
  * [71] C. Voivret, F. Radja??, J.-Y. Delenne, M. S. E. Youssoufi, Space-filling properties of polydisperse granular media, Physical Review E : Statistical, Nonlinear, and Soft Matter Physic 76 (2007) 021301.
  * [72] C. Voivret, F. Radjai, J.-Y. Delenne, M. S. E. Youssoufi, Multiscale force networks in highly polydisperse granular media, Physical Review Letter 102 (2009) 178001.
  * [73] P. Y. Liu, R. Y. Yang, A. B. Yu, The effect of liquids on radial segregation of granular mixtures in rotating drums, Granular Matter 15 (4) (2013) 427–436.


# Providing Previously Unseen Users Fair Recommendations Using Variational
Autoencoders

Bjørnar Vassøy, Helge Langseth and Benjamin Kille
Norwegian University of Science and Technology

###### Abstract

An emerging definition of fairness in machine learning requires that models
are oblivious to demographic user information, e.g., a user’s gender or age
should not influence the model. Personalized recommender systems are
particularly prone to violating this definition through their explicit user
focus and user modelling. Explicit user modelling is also an aspect that makes
many recommender systems incapable of providing hitherto unseen users with
recommendations. We propose novel approaches for mitigating discrimination in
Variational Autoencoder-based recommender systems by limiting the encoding of
demographic information. The approaches are capable of, and evaluated on,
providing users that are not represented in the training data with fair
recommendations.

##  1 Introduction

Fairness in recommender systems is becoming a popular and diverse research
field. Burke (2017) formalized the multi-stakeholder nature of the
recommendation setting: Producers are interested in mitigating popularity bias
(Ahanger et al., 2022) to ensure their products are given the exposure they
deserve. Consumers expect to be treated similarly regardless of demographic
attributes like age, race, and gender. In addition to the stakeholder
perspectives, there is no single definition of what constitutes a fair
recommendation. A consumer-side perspective focusing on discrimination of
demographic user groups may require that similar ratings are estimated for
each group (Kamishima et al., 2018), that each group receives similar
recommendations (Farnadi et al., 2018), that each group are equally satisfied
with their recommendations (Yao and Huang, 2017), or that model
representations do not correlate with the user groups. This research focuses
on the latter notion (Neutral Representations) while also evaluating whether
similar recommendations are given to different demographic groups
(Recommendation Parity).

A mostly unexplored perspective of fair recommender systems relates to
recommending for users that were not part of the training data. A user’s first
impression of a platform may decide if they will continue using it. A man who
aspires to be a florist and likes romance movies may feel stereotypes and be
discouraged from further interaction if he is first recommended physical
labour careers or action movies. Unlike contemporary fair recommender systems,
our research focuses on introducing fairness in a model architecture that can
recommend for all users, including users not represented in training data. We
propose Variational Autoencoder (VAE) approaches that only require a list of
items a user has interacted with and no pre-trained user representations. The
approaches may replace contextual and item-to-item recommender systems used to
onboard new users or serve as a complete recommendation platform that does not
require frequent model updates.

One goal of this research is to provide more insight into the competitiveness
of VAE-based recommender systems given their limited use in recommender system
research. Given this, we wish to verify if our VAE-based models can fairly
process unseen users and if the same models produce state-of-the-art fair
recommendations for said users. Formalized research questions are as follows:
1. Are VAE-based recommender systems competitive? and 2. Can the encoded
demographic information in VAE latent states be reduced when processing new
users?

##  2 Related Work

Multiple methods have been proposed for filtering out sensitive information
embedded in model representations and parameters. Many of these methods train
adversarial models tasked with classifying sensitive attributes given
representations belonging to the recommender systems. The recommender systems
are then penalized for encoding sensitive information by adding additional
objectives of fooling the adversarial models. This strategy has been applied
directly to latent factors of factorization models (Resheff et al., 2019; Xu
et al., 2021; Wu et al., 2021a), indirectly to train attribute filters (Bose
and Hamilton, 2019), and for filtering graph embeddings (Wu et al., 2021b; Liu
et al., 2022). Another strategy encourages representation neutrality by
optimizing for representations that are orthogonal to sensitive dimensions in
the representation space (Wu et al., 2021a; Islam et al., 2021). There are
also examples of introducing neutrality by adjusting sampling schemes applied
while training representations (Rahman et al., 2019; Li et al., 2022), and
methods for causally isolating sensitive information to specific model factors
that are replaced or dropped during inference (Buyl and Bie, 2020; Frisch et
al., 2021).

The framework proposed by Li et al. (2023) shares the most high-level
similarities with this work in considering Representation Neutrality fairness
optimized through adversarial methods and in being capable of recommending for
unseen users. The latter is achieved by training an auxiliary mapping function
to map new users to the representation space of the trained model.

Creager et al. (2019) propose a VAE-based classification model that allows its
users to specify which sensitive attributes it should be oblivious to. One of
the key contributions is to isolate all sensitive user information in one part
of the latent representation. The sensitive part of the latent representation
is subjected to a secondary task of classifying the user’s sensitive
attributes, and an approximated KL-divergence measure is minimized to
encourage independence between the two parts of the latent representation.

##  3 Background

###  3.1 Variational Autoencoders

The Variational Autoencoder (VAE) is a variational Bayesian model initially
proposed by Kingma and Welling (2014) that has since seen wide application
within representation learning and image generation. The vanilla VAE is posed
as a graphical model where the observed data $\boldsymbol{x}$ depends on a
latent variable $\boldsymbol{z}$. The final objective of maximizing the
Evidence Lower Bound (ELBO) can be derived from minimizing the KL-divergence
from the variational distribution $q(\boldsymbol{z}|\boldsymbol{x})$ to the
true posterior $p(\boldsymbol{z}|\boldsymbol{x})$.

$\text{ELBO}(q,p)=\mathbb{E}_{q(\boldsymbol{z}|\boldsymbol{x})}\left[\log
p(\boldsymbol{x}|\boldsymbol{z})\right]-D_{KL}\left[q(\boldsymbol{z}|\boldsymbol{x})||p(\boldsymbol{z})\right],$
(1)

where $p(\boldsymbol{x}|\boldsymbol{z})$ is the likelihood and
$p(\boldsymbol{z})$ is the latent prior. Neural networks typically
parameterize $p(\boldsymbol{x}|\boldsymbol{z})$ and
$q(\boldsymbol{z}|\boldsymbol{x})$.

Higgins et al. (2017) proposed the $\beta$-VAE which adds a $\beta$ factor
associated with the KL-divergence term. In particular, they explore $\beta>1$
to produce disentangled encodings, i.e., the increased focus on independence
between dimensions of $\boldsymbol{z}$ leads to semantically different
concepts encoded in them.

Similarly, Kim and Mnih (2018) adds an optimization term for increased
disentanglement that penalizes dependency between dimensions of the latent
representation.

$\text{Loss}_{\text{FactorVAE}}(q,p)=\text{ELBO}(q,p)-\gamma
D_{KL}\left[q(\boldsymbol{z})||\prod_{i}q(\boldsymbol{z}_{i})\right]$ (2)

###  3.2 VAE-based Recommender System

Liang et al. (2018) propose applying a VAE as a Recommender System by having
$\boldsymbol{x}$ represent user interaction history, e.g., items that the user
has interacted with or rated. Recommendations are extracted from the fuzzy
reconstruction of the user interaction history based on values assigned to
items the user has not previously interacted with.

##  4 Methodology

The key idea behind our model is to leverage the bottleneck characteristics of
the VAE to encourage the reduction of sensitive user information such that the
provided recommendations are minimally influenced by such information. The
probabilistic nature and low dimensionality of the latent representation can
simplify the objective of filtering out sensitive information. The main
proposed model setups are illustrated in Figure 1.

Figure 1: Illustration of the Split Latent model setups. The encoder is
dynamically designed and is, in practice, implemented as two separate
encoders, one for $\boldsymbol{z}$ and one for $\boldsymbol{b}$. Key details
are that no explicit sensitive information $\boldsymbol{s}$ is provided or
required during inference, and the sensitive part of the latent representation
$\boldsymbol{b}$ is not used for recommendation.

###  4.1 Base Model

The underlying recommender system is based on the model proposed by Liang et
al. (2018). In addition to posing the recommendation task as the objective of
reconstructing the users’ interaction histories, Liang et al. (2018) propose
multiple extensions. They apply dropout on the input during training (noisy
VAE) to improve generalization. They explore Gaussian and Logistic likelihood
$p(\boldsymbol{x}|\boldsymbol{z})$ before settling on the Multinomial for its
performance and nice properties. The multinomial likelihood does not
explicitly penalize probability density allocated to the items the user has
not interacted with, which in turn avoids the assumption of many other options
in that these are all items that the user dislike and that should be allocated
zero probability.

$\log
p(\boldsymbol{x}|\boldsymbol{z})=\log\left[\prod_{i}p(\boldsymbol{x}_{i}|\boldsymbol{z})^{\boldsymbol{x}_{i}}\right]=\sum_{i}\boldsymbol{x}_{i}\text{log}p(\boldsymbol{x}_{i}|\boldsymbol{z}),\text{where}\sum_{i}p(\boldsymbol{x}_{i}|\boldsymbol{z})=1$
(3)

A specific item $\boldsymbol{x}_{i}$’s contribution is zero if it is not found
in the user’s interaction history, regardless of what the decoder
parameterizes the item’s probability $p(\boldsymbol{x}_{i}|\boldsymbol{z})$ to
be. Unlike Higgins et al. (2017), Liang et al. (2018) explore $\beta<1$,
citing that generation has limited applications in the recommender system
setting, and identify $\beta=0.2$ as a good candidate.

The model of Liang et al. (2018) was altered slightly in this work: The
proposed $\beta$-annealing strategy was dropped since it did not yield
noticeable improvements, the Hyperbolic Tangent activation was switched out
with the SELU (Klambauer et al., 2017), and better results were achieved when
reducing the dimensionality of the latent state. The latent dimension was set
to 64 for the baseline setup and 24 for all fairness setups, as opposed to 200
in Liang et al. (2018). Further, $\beta$ was set to 1 for the extended setups
as it was found to synergize better with the fairness extensions.

###  4.2 Adversarial Setup

The adversarial setup is considered a baseline setup and couples the base
model with an adversarial model tasked with classifying sensitive attributes
from the latent state $\boldsymbol{z}$. This is one of the extensions explored
by Borges and Stefanidis (2022). Insight from the adversarial model help the
main model avoid encoding sensitive information.

###  4.3 Split Latent Setups

The split latent state setups bisect the latent state into one part for
encoding sensitive user information $\boldsymbol{b}$ and another part that is
free for sensitive information $\boldsymbol{z}$. Fair recommendations are
produced by decoding $\boldsymbol{z}$, while $\boldsymbol{b}$ is discarded
during evaluation. The motivation behind the bisected latent state is to
leverage the sensitive information in $\boldsymbol{b}$ to inform the model of
the information that should not be encoded in $\boldsymbol{z}$. An
inspirational split latent setup is proposed by Creager et al. (2019). Our
approaches differ in that we consider recommendation rather than
classification, we do not limit the encoding of each sensitive attribute to
single dimensions in $\boldsymbol{b}$, and we posit isotropic Gaussian priors
to $\boldsymbol{b}$.

A classification task is coupled with $\boldsymbol{b}$ through another
decoder, which will be referred to as the Sensitive Decoder, with the goal of
inferring the user’s sensitive attributes. In both Split Latent setups, binary
sensitive attributes were considered, and the re-classification was optimized
using cross-entropy. The choice of an isotropic Gaussian prior
$p(\boldsymbol{z},\boldsymbol{b})$ will inherently optimize the VAE to produce
independent dimensions of $[\boldsymbol{z}\;\boldsymbol{b}]$, but this is
supplemented with an explicit term for penalizing correlation between
$\boldsymbol{z}$ and $\boldsymbol{b}$. The full objective to be maximized is

$\begin{split}\text{Split}&\text{LatentObj}(q,p)=\\\
&\mathbb{E}_{q(\boldsymbol{z},\boldsymbol{b}|\boldsymbol{x})}\left[\log
p(\boldsymbol{x}|\boldsymbol{z})+\alpha\log
p(\boldsymbol{s}|\boldsymbol{b})\right]-\beta
D_{KL}\left[q(\boldsymbol{z},\boldsymbol{b}|\boldsymbol{x})||p(\boldsymbol{z},\boldsymbol{b})\right]-\gamma
D_{KL}\left[q(\boldsymbol{z},\boldsymbol{b})||q(\boldsymbol{z})q(\boldsymbol{b})\right],\end{split}$
(4)

where $\beta$, $\alpha$ and $\gamma$ are hyperparameters for adjusting the
influence of terms,
$\mathbb{E}_{q(\boldsymbol{z}|\boldsymbol{x})}\left[\alpha\text{log}p(\boldsymbol{s}|\boldsymbol{b})\right]$
comprise the re-classification objective of sensitive attribute
$\boldsymbol{s}$ given $\boldsymbol{b}$, and
$D_{KL}\left[q(\boldsymbol{z},\boldsymbol{b})||q(\boldsymbol{z})q(\boldsymbol{b})\right]$
is the term introduced for further penalizing correlation between
$\boldsymbol{z}$ and $\boldsymbol{b}$. The latter is defined and implemented
in two different ways.

Generative Adversarial Network (GAN) KL. Inspiration has been taken from Kim
and Mnih (2018), who proposed approximating the KL divergence from the
aggregate posterior
$\hat{q}(\boldsymbol{z})=\mathbb{E}_{\boldsymbol{x}\in\mathcal{D}_{s}}\left[q(\boldsymbol{z}|\boldsymbol{x})\right]$,
where $\mathcal{D}_{s}$ is the sampled minibatch, to a factorization over each
dimension using an adversarial model. This measure is approximated by an
adversarial model trained to tell the original latent representations apart
from factorized ones where dimensions have been shuffled across the minibatch.
The equivalent for $q(\boldsymbol{z})q(\boldsymbol{b})$ is to shuffle
$\boldsymbol{b}$ in its entirety across the minibatch.

Empiric KL. GAN approaches are often hard to optimize due to the moving
objectives and the balancing of two competing models. An alternative approach
is based on the analytic KL-divergence of Gaussian distributions. The empiric
covariance of a minibatch of latent representations replaces the covariances
of $q(\boldsymbol{z},\boldsymbol{b})$, while the covariances of
$q(\boldsymbol{z})q(\boldsymbol{b})$ are replaced by the same empiric
covariance matrix where all covariances between $\boldsymbol{z}$ and
$\boldsymbol{b}$ are set to zero.

$D_{KL}\left[q(\boldsymbol{z},\boldsymbol{b})||q(\boldsymbol{z})q(\boldsymbol{b})\right]=\frac{1}{2}\left[\log\frac{|\hat{\Sigma}_{2}|}{|\hat{\Sigma}_{1}|}-d+\text{tr}\left\\{\hat{\Sigma}_{2}^{-1}\hat{\Sigma}_{1}\right\\}\right],\text{when
}\hat{\mu}_{1}=\hat{\mu}_{2},$ (5)

where $\hat{\mu}$ are empiric means, $\hat{\Sigma}$ are empiric covariances,
$d$ is the number of dimensions in the latent state, and $\text{tr}\\{\\}$ is
the trace operator. The covariance matrix $\hat{\Sigma}_{2}$ is block-
diagonal, so its inverse can be computed block by block.

##  5 Results

###  5.1 Experimental Setup

Two established datasets were used to conduct the experiments. Unlike
contemporary work, both datasets were split into training, validation and test
datasets containing disjoint sets of users to accommodate the setting of
providing unseen users with recommendations. Common for both datasets was the
choice of two binary sensitive attributes: age and gender. The datasets
provide gender as ‘male’, ‘female’, or various forms for
‘missing’/‘undefined’/‘other’, but the application of labels other than male
and female was deemed too inconsistent to make out one or more additional
sensitive labels. Age was also made binary to complement the gender attribute,
and a threshold age of 35 was chosen. Users that miss either sensitive
attribute are filtered out since these are required during training and
evaluation.
MovieLens1M: The 1 million version of the MovieLens dataset (Harper and
Konstan, 2015) is the most applied dataset for evaluating consumer-side
fairness in recommender systems. The raw dataset contains 1 million movie
ratings and is the largest MovieLens dataset that provides user attributes
that are considered sensitive. Following established practices for converting
ratings into implicit feedback data, ratings of 4 or 5 out of 5 were labelled
1 while lower ratings were set to 0.
LastFM 2B: The LFM-2b dataset (Melchiorre et al., 2021) comprises 2 billion
listening events collected on the lastFM platform. This dataset was processed
to focus on the most recent and relevant data, as well as reducing the number
of items. The most recent two years of data were extracted, and the objective
was set to recommend artists rather than albums or songs. Finally, artists
with fewer than 110 listening events were filtered out. Key statistics of the
datasets are summarized in Table 1.

Table 1: Key dataset statistics. Dataset | # Users | #Female/#Male | #Senior/#Young | # Items | # Records
---|---|---|---|---|---
MovieLens | 6k |  1.7k/4.3k |  1.4k/4.6k |  3.5k |  575k
LastFM |  14.6k |  3k/11.5k |  1.3k/13.3k |  22k |  8750k

####  5.1.1 Metrics

NDCG@k: NDCG@k was used as the recommendation utility metric. NDCG is a
popular ranking metric in recommender systems that rewards ranking relevant
items high. For all experiments, k was set to 10.
AUC: AUC(Area Under the ROC Curve) is a metric commonly used for
classification that considers the True Positive Rate and the False Positive
Rate for all possible threshold values. For fair recommender systems it has
been used to evaluate how well sensitive information has been filtered out of
representations, which is its role in this work. For this particular objective
and binary sensitive attributes, the perfect AUC score is 0.5 which indicates
that all sensitive information has been filtered out. AUC is measured by
training an auxiliary classification model on the model representations. One
challenge with this metric is that it cannot be applied to the SLIM baseline.
$\boldsymbol{\chi^{2}}$-statistic: $\chi^{2}$-test can be used to estimate the
probability that two independent samples were drawn from the same
distribution. For the considered setting, it is natural to compare how often
items occur in recommendations given to different sensitive groups, e.g.,
young and senior users. Top recommendations are not independent, rendering the
test unreliable, and we instead focus on the statistic. The top 100
recommendations given to each user were selected to aggregate the contingency
tables of each sensitive attribute. The number of items considered was set for
each dataset such that each cell had an expected value of at least 3 in all
setups since the long-tail nature of recommendation yields a lot of cells with
low expectations and since the statistic is highly sensitive to such cells.
Expectations are adjusted according to the number of users who can be
recommended each item, i.e., users who have not already interacted with it.
One issue with considering the top recommendations observations is that we
lose the ordering information in the ranking, where rank 1 is assigned more
confidence than rank 100.
Kendall-Tau distance: Kendall-Tau distance is a distance metric between two
ordered lists. The original distance is undefined when items only occur in one
of the lists, so an extension designed for the recommender
setting111https://godatadriven.com/blog/using-kendalls-tau-to-compare-
recommendations/ is applied. The extension is shown to have some intuitive
properties and will output values from 1, being a perfect match, to -1, when
the original recommendation lists contain disjointed sets of items. This
metric was also considered for top-$k$ recommendations given to sensitive
groups, where $k$ is set to 100. The recommendations assigned to each
sensitive group were aggregated over individual user recommendations while
applying a rank discounting scheme to give higher importance to highly ranked
items and to serve as a shared normalization strategy for the different
models.

####  5.1.2 Models

Since no implementations of viable fair baselines have been identified,
comparisons are made with non-fair baselines and focus on how the different
proposed fairness extensions of a base VAE recommender system impact
recommendation utility and fairness measures. All code is publically available
on GitHub222https://github.com/BjornarVass/fair-vae-rec.
SLIM: SLIM (Ning and Karypis, 2011) is an established baseline recommender
system that culminates in an item-to-item parameter matrix. The item-to-item
nature allows it to recommend for users not considered during model fitting.
VAErec: The model proposed by Liang et al. (2018) with minor alterations as
specified in Section 4.1. This model is not trained with any fairness
objectives and serves as the base model that is extended to make out the other
VAE models.
VAEadv: VAERec extended with an adversarial model that filters sensitive
information from the latent representation.
VAEgan: VAERec extended with bisected latent representation for isolating and
filtering out sensitive information. Independence between sections is
optimized using an adversarial model that approximates a KL-divergence term.
VAEemp: VAERec extended with bisected latent representation for isolating and
filtering out sensitive information. Independence between sections is
optimized using an analytic KL-divergence term with empiric covariances.

###  5.2 Main Results

VAErec achieves comparable NDCG as SLIM on the MovieLens dataset, but SLIM
performs better on the LastFM dataset. This is contrary to the results of
Liang et al. (2018) which showed improved NDCG@100 over SLIM on two movie
datasets, of which the results on the 20 million version of MovieLens were
successfully reproduced using our implementation. It is unclear which factors
affect the performance of the VAE-recommender, particularly since our changes
further improved the original model and considering that two datasets from the
same source conflict in performance.

All fairness extensions are shown to improve the AUC scores significantly. In
particular, the AUC scores of VAErec on the MovieLens dataset indicate that
the latent aptly encodes user gender and age with AUCs above 0.8, whereas
VAEemp reduces this to 0.65-0.63 at the cost of 0.035 NDCG ($\approx 11\%$).
VAEemp outperforms the other options when considering NDCG and Representation
Neutrality and manages to filter out age better than gender despite the base
model inherently encoding more age information. The secondary fairness metrics
suggest that all fair extensions produce more similar recommendations for
users of different sensitive groups than SLIM and VAErec. The
$\chi^{2}$-statistic of VAErec is worse than that of SLIM, but the fair
extensions all significantly outperform SLIM. Interestingly, VAEemp is the
worst-performing fair extension on these metrics, which suggests they do not
perfectly reflect the AUC metric.

Table 2: MovieLens results. Model | NDCG@10$\uparrow$ | AUC G$\downarrow$ | AUC A$\downarrow$ |  $\boldsymbol{\chi^{2}}$@100 G$\downarrow$ |  $\boldsymbol{\chi^{2}}$@100 A$\downarrow$ | K.T@100 G$\uparrow$ | K.T@100 A$\uparrow$
---|---|---|---|---|---|---|---
SLIM |  0.328$\pm$0.009 | - | - |  2285$\pm$280.1 |  2198$\pm$237.6 |  0.476$\pm$0.075 |  0.448$\pm$0.045
VAErec |  0.321$\pm$0.008 |  0.804$\pm$0.024 |  0.859$\pm$0.019 |  2990$\pm$415.9 |  2636$\pm$359.3 |  0.559$\pm$0.054 |  0.537$\pm$0.035
VAEadv |  0.280$\pm$0.008 |  0.678$\pm$0.036 |  0.675$\pm$0.043 |  1121$\pm$273.4 |  904.6$\pm$194.5 |  0.820$\pm$0.025 |  0.792$\pm$0.038
VAEgan |  0.277$\pm$0.010 |  0.687$\pm$0.037 |  0.695$\pm$0.050 |  1054$\pm$232.6 |  852.5$\pm$208.0 |  0.854$\pm$0.036 |  0.841$\pm$0.029
VAEemp |  0.286$\pm$0.008 |  0.652$\pm$0.032 |  0.629$\pm$0.041 |  1355$\pm$302.7 |  1151$\pm$228.8 |  0.804$\pm$0.033 |  0.770$\pm$0.043

The LastFM dataset offers different insights on minority groups since the
sensitive groups are very skewed, e.g., only 10% of the users are considered
senior. This seems to be reflected in a smaller improvement of Representation
Neutrality achieved by the fair extensions compared with their performance on
the MovieLens dataset. The improved gender AUC is comparable, but the base AUC
of VAErec is lower. All three fair extensions performed similarly on gender,
but VAEemp struggles more than the other two on age. The extensions’ lesser
improvements of AUC over VAErec come at a smaller relative reduction in NDCG
than the one seen on the Movielens dataset (¡ 4% for VAEadv and VAEemp).

Kendall Tau on LastFM is significantly better for all VAE-based models, and
the fair models all outperform VAErec. On the other hand, only VAEadv achieve
a $\chi^{2}$-statistic that is better than that of SLIM for this dataset. One
confounding factor is that SLIM is observed to provide more diverse
recommendations, with roughly 60% of recommendations being the top 10% popular
items, vs roughly 85% for VAEemp, meaning that the 750 most popular items
considered in the $\chi^{2}$-statistic of LastFM cover far fewer of the total
SLIM recommendations. When comparing the VAE-based models, it is clear that
the fair extensions succeed in reducing the very large initial
$\chi^{2}$-statistics of VAErec.

Table 3: LastFM results. Model | NDCG@10$\uparrow$ | AUC G$\downarrow$ | AUC A$\downarrow$ |  $\boldsymbol{\chi^{2}}$@100 G$\downarrow$ |  $\boldsymbol{\chi^{2}}$@100 A$\downarrow$ | K.T@100 G$\uparrow$ | K.T@100 A$\uparrow$
---|---|---|---|---|---|---|---
SLIM |  0.572$\pm$0.010 | - | - |  2681$\pm$481.7 |  2164$\pm$292.6 |  0.238$\pm$0.099 |  -0.032$\pm$0.061
VAErec |  0.512$\pm$0.010 |  0.759$\pm$0.014 |  0.821$\pm$0.023 |  3693$\pm$661.9 |  2840$\pm$485.3 |  0.608$\pm$0.046 |  0.440$\pm$0.073
VAEadv |  0.492$\pm$0.010 |  0.666$\pm$0.019 |  0.731$\pm$0.025 |  2441$\pm$546.6 |  2068$\pm$416.8 |  0.757$\pm$0.041 |  0.655$\pm$0.043
VAEgan |  0.484$\pm$0.014 |  0.674$\pm$0.035 |  0.738$\pm$0.031 |  2596$\pm$656.8 |  2191$\pm$601.3 |  0.775$\pm$0.054 |  0.683$\pm$0.061
VAEemp |  0.495$\pm$0.010 |  0.663$\pm$0.022 |  0.757$\pm$0.027 |  2924$\pm$577.4 |  2528$\pm$481.4 |  0.725$\pm$0.045 |  0.588$\pm$0.046

###  5.3 Sampling Feature

VAE-based models typically only consider the parameterized mean of the
variational distribution during inference, i.e., it is deterministic. In a
setting where fairness is of utmost importance, one can leverage the
parameterized mean and variation to sample latent states that are inherently
noisy. Sampled latent states are typically fairer, i.e., more neutral, but
produce less accurate recommendations. Thus, VAE-based recommender systems can
dynamically offer two different modes depending on how the user values the
performance and fairness tradeoff.

The sampling feature is compared for different values of the hyperparameter
$\beta$ since it directly controls the loss term regulates the parameterized
distribution. The default setting of $\beta=1.0$ yields the biggest difference
in results when applying sampled or deterministic latent states. Sampling with
$\beta=1.0$ resulted in the best fairness scores but also the worst NDCG
score. On the other hand, $\beta$ set to 0.2 and 0.6 produced marginally
better NDCG and AUC for the deterministic mode with $\beta=0.2$ coming out on
top. Reducing $\beta$ appears to improve the fairness of the deterministic
mode at the loss of fairness in the sampling mode. Notably, sampling with
$\beta=0.2$ and not sampling with $\beta=1$ achieves similar NDCG, but the
AUCs achieved by the former is noticably lower. This suggests that coupling
the choice of $\beta$ and the sampling strategy may yield good settings for
scenarios where one metric is assigned a strict upper or lower bound
constraint. Large $\beta$ may be ideal in dynamic scenarios where users can
choose to turn on sampling to improve fairness.

Table 4: Results from different $\beta$ settings, with and without sampling latent representations. Model | NDCG@10$\uparrow$ | AUC G$\downarrow$ | AUC A$\downarrow$
---|---|---|---
VAEemp $\beta=1.0$ |  0.286$\pm$0.008 |  0.651$\pm$0.032 |  0.629$\pm$0.041
sampled |  0.256$\pm$0.009 |  0.595$\pm$0.035 |  0.562$\pm$0.026
VAEemp $\beta=0.6$ |  0.292$\pm$0.009 |  0.652$\pm$0.029 |  0.615$\pm$0.042
sampled |  0.269$\pm$0.007 |  0.603$\pm$0.036 |  0.573$\pm$0.038
VAEemp $\beta=0.2$ |  0.292$\pm$0.008 |  0.640$\pm$0.044 |  0.607$\pm$0.030
sampled |  0.279$\pm$0.008 |  0.619$\pm$0.042 |  0.587$\pm$0.038

##  6 Conclusion and Future Work

Relating to the first research question, this research indicates that VAE-
based recommenders consistently perform well but are not universally
competitive. While few architectures excel in all scenarios, this insight may
motivate research into other model types capable of providing new users with
fair recommendations. To answer the second research question, we have shown
that these models can be successfully applied and extended to significantly
limit the demographic information encoded in the latent state, meaning that
they can provide new and established users alike with fairer recommendations.
The improvement in fairness comes with a minor deterioration of recommendation
utility, as is seen in similar research. The VAE also offers a means for
further obfuscating user representation through parameterized latent state
variance. For the primary fairness definition of Neutral Representation, the
proposed extensions significantly outperform the base model. The models also
perform well on the secondary fairness definition of Recommendation Parity
where the proposed VAE extensions either outperform or match relevant models
on two metrics.

For future work, it could be interesting to evaluate the model on other
datasets representing different recommendation settings. It would also be
desirable to compare with other models that optimize for fairness when
recommending for new users. Such models could be based on contextual
recommender systems or other architectures that do not learn explicit user
representations.

## Acknowledgments

This publication has been partly funded by the SFI NorwAI, (Centre for
Research-based Innovation, 309834). The authors gratefully acknowledge the
financial support from the Research Council of Norway and the partners of the
SFI NorwAI.

## References

  * Ahanger et al. [2022] Abdul Basit Ahanger, Syed Wajid Aalam, Muzafar Rasool Bhat, and Assif Assad.  Popularity bias in recommender systems - a review.  In Valentina E. Balas, G. R. Sinha, Basant Agarwal, Tarun Kumar Sharma, Pankaj Dadheech, and Mehul Mahrishi, editors, _Emerging Technologies in Computer Engineering: Cognitive Computing and Intelligent IoT_ , pages 431–444, Cham, 2022. Springer International Publishing.  ISBN 978-3-031-07012-9.
  * Borges and Stefanidis [2022] Rodrigo Borges and Kostas Stefanidis.  F2VAE: a framework for mitigating user unfairness in recommendation systems.  In _Proceedings of the 37th ACM/SIGAPP Symposium on Applied Computing_ , pages 1391–1398, Virtual Event, April 2022. ACM.  ISBN 978-1-4503-8713-2.  doi: 10.1145/3477314.3507152.
  * Bose and Hamilton [2019] Avishek Bose and William Hamilton.  Compositional Fairness Constraints for Graph Embeddings.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 715–724, Long Beach Convention Center, Long Beach, June 2019\. PMLR.  URL https://proceedings.mlr.press/v97/bose19a.html.
  * Burke [2017] Robin Burke.  Multisided fairness for recommendation.  _CoRR_ , abs/1707.00093, 2017.  URL http://arxiv.org/abs/1707.00093.
  * Buyl and Bie [2020] Maarten Buyl and Tijl De Bie.  DeBayes: a Bayesian Method for Debiasing Network Embeddings.  In _Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 1220–1229, Virtual, 2020. PMLR.  URL http://proceedings.mlr.press/v119/buyl20a.html.
  * Creager et al. [2019] Elliot Creager, David Madras, Joern-Henrik Jacobsen, Marissa Weis, Kevin Swersky, Toniann Pitassi, and Richard Zemel.  Flexibly fair representation learning by disentanglement.  In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pages 1436–1445, Long Beach Convention Center, Long Beach, 09–15 Jun 2019. PMLR.  URL https://proceedings.mlr.press/v97/creager19a.html.
  * Farnadi et al. [2018] Golnoosh Farnadi, Pigi Kouki, Spencer K. Thompson, Sriram Srinivasan, and Lise Getoor.  A fairness-aware hybrid recommender system.  _CoRR_ , abs/1809.09030, 2018.  URL http://arxiv.org/abs/1809.09030.
  * Frisch et al. [2021] Gabriel Frisch, Jean-Benoist Leger, and Yves Grandvalet.  Co-clustering for Fair Recommendation.  In Michael Kamp, Irena Koprinska, Adrien Bibal, Tassadit Bouadi, and Benoît Frénay, editors, _Machine Learning and Principles and Practice of Knowledge Discovery in Databases_ , pages 607–630, Cham, 2021\. Springer International Publishing.  ISBN 978-3-030-93736-2.
  * Harper and Konstan [2015] F. Maxwell Harper and Joseph A. Konstan.  The movielens datasets: History and context.  _ACM Trans. Interact. Intell. Syst._ , 5(4), dec 2015.  ISSN 2160-6455.  doi: 10.1145/2827872.  URL https://doi.org/10.1145/2827872.
  * Higgins et al. [2017] Irina Higgins, Loïc Matthey, Arka Pal, Christopher P. Burgess, Xavier Glorot, Matthew M. Botvinick, Shakir Mohamed, and Alexander Lerchner.  beta-vae: Learning basic visual concepts with a constrained variational framework.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_ , Toulon, France, 2017. OpenReview.net.  URL https://openreview.net/forum?id=Sy2fzU9gl.
  * Islam et al. [2021] Rashidul Islam, Kamrun Naher Keya, Ziqian Zeng, Shimei Pan, and James Foulds.  Debiasing Career Recommendations with Neural Fair Collaborative Filtering.  In _Proceedings of the Web Conference 2021_ , pages 3779–3790, Ljubljana Slovenia, April 2021. ACM.  ISBN 978-1-4503-8312-7.  doi: 10.1145/3442381.3449904.
  * Kamishima et al. [2018] Toshihiro Kamishima, Shotaro Akaho, Hideki Asoh, and Jun Sakuma.  Recommendation Independence.  In Sorelle A. Friedler and Christo Wilson, editors, _Conference on Fairness, Accountability and Transparency, FAT 2018, 23-24 February 2018, New York, NY, USA_ , volume 81 of _Proceedings of Machine Learning Research_ , pages 187–201, New York, USA, 2018. PMLR.  URL http://proceedings.mlr.press/v81/kamishima18a.html.
  * Kim and Mnih [2018] Hyunjik Kim and Andriy Mnih.  Disentangling by factorising.  In Jennifer Dy and Andreas Krause, editors, _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pages 2649–2658, Stockholm, Sweden, 10–15 Jul 2018. PMLR.  URL https://proceedings.mlr.press/v80/kim18b.html.
  * Kingma and Welling [2014] Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  In Yoshua Bengio and Yann LeCun, editors, _2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings_ , Banff, AB, Canada, 2014. ICLR.  URL http://arxiv.org/abs/1312.6114.
  * Klambauer et al. [2017] Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter.  Self-normalizing neural networks.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30, Montreal, Canada, 2017. Curran Associates, Inc.  URL https://proceedings.neurips.cc/paper_files/paper/2017/file/5d44ee6f2c3f71b73125876103c8f6c4-Paper.pdf.
  * Li et al. [2022] Cheng-Te Li, Cheng Hsu, and Yang Zhang.  Fairsr: Fairness-aware sequential recommendation through multi-task learning with preference graph embeddings.  _ACM Trans. Intell. Syst. Technol._ , 13(1), feb 2022.  ISSN 2157-6904.  doi: 10.1145/3495163.
  * Li et al. [2023] Yunqi Li, Dingxian Wang, Hanxiong Chen, and Yongfeng Zhang.  Transferable fairness for cold-start recommendation, 2023.
  * Liang et al. [2018] Dawen Liang, Rahul G. Krishnan, Matthew D. Hoffman, and Tony Jebara.  Variational autoencoders for collaborative filtering.  In _Proceedings of the 2018 World Wide Web Conference_ , WWW ’18, page 689–698, Republic and Canton of Geneva, CHE, 2018. International World Wide Web Conferences Steering Committee.  ISBN 9781450356398.  doi: 10.1145/3178876.3186150.  URL https://doi.org/10.1145/3178876.3186150.
  * Liu et al. [2022] Haifeng Liu, Hongfei Lin, Wenqi Fan, Yuqi Ren, Bo Xu, Xiaokun Zhang, Dongzhen Wen, Nan Zhao, Yuan Lin, and Liang Yang.  Self-supervised learning for fair recommender systems.  _Applied Soft Computing_ , 125:109126, August 2022.  ISSN 1568-4946.  doi: 10.1016/j.asoc.2022.109126.
  * Melchiorre et al. [2021] Alessandro B. Melchiorre, Navid Rekabsaz, Emilia Parada-Cabaleiro, Stefan Brandl, Oleg Lesota, and Markus Schedl.  Investigating gender fairness of recommendation algorithms in the music domain.  _Information Processing & Management_, 58(5):102666, 2021.  ISSN 0306-4573.  doi: https://doi.org/10.1016/j.ipm.2021.102666.  URL https://www.sciencedirect.com/science/article/pii/S0306457321001540.
  * Ning and Karypis [2011] Xia Ning and George Karypis.  SLIM: sparse linear methods for top-n recommender systems.  In Diane J. Cook, Jian Pei, Wei Wang, Osmar R. Zaïane, and Xindong Wu, editors, _11th IEEE International Conference on Data Mining, ICDM 2011, Vancouver, BC, Canada, December 11-14, 2011_ , pages 497–506, Vancouver, Canada, 2011. IEEE Computer Society.  doi: 10.1109/ICDM.2011.134.  URL https://doi.org/10.1109/ICDM.2011.134.
  * Rahman et al. [2019] Tahleen Rahman, Bartlomiej Surma, Michael Backes, and Yang Zhang.  Fairwalk: Towards Fair Graph Embedding.  In _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence_ , pages 3289–3295, Macao, China, August 2019. International Joint Conferences on Artificial Intelligence Organization.  ISBN 978-0-9992411-4-1.  doi: 10.24963/ijcai.2019/456.
  * Resheff et al. [2019] Yehezkel S. Resheff, Yanai Elazar, Moni Shahar, and Oren Sar Shalom.  Privacy and Fairness in Recommender Systems via Adversarial Training of User Representations.  In Maria De Marsico, Gabriella Sanniti di Baja, and Ana L. N. Fred, editors, _Proceedings of the 8th International Conference on Pattern Recognition Applications and Methods, ICPRAM 2019, Prague, Czech Republic, February 19-21, 2019_ , pages 476–482, Prague, Czech Republic, 2019. SciTePress.  doi: 10.5220/0007361204760482.
  * Wu et al. [2021a] Chuhan Wu, Fangzhao Wu, Xiting Wang, Yongfeng Huang, and Xing Xie.  Fairness-aware news recommendation with decomposed adversarial learning.  In _Thirty-Fifth AAAI Conference on Artificial Intelligence, AAAI 2021, Thirty-Third Conference on Innovative Applications of Artificial Intelligence, IAAI 2021, The Eleventh Symposium on Educational Advances in Artificial Intelligence, EAAI 2021, Virtual Event, February 2-9, 2021_ , pages 4462–4469, Virtual, 2021a. AAAI Press.  URL https://ojs.aaai.org/index.php/AAAI/article/view/16573.
  * Wu et al. [2021b] Le Wu, Lei Chen, Pengyang Shao, Richang Hong, Xiting Wang, and Meng Wang.  Learning Fair Representations for Recommendation: A Graph-Based Perspective.  In _Proceedings of the Web Conference 2021_ , WWW ’21, pages 2198–2208, New York, NY, USA, 2021b. Association for Computing Machinery.  ISBN 978-1-4503-8312-7.  doi: 10.1145/3442381.3450015.  event-place: Ljubljana, Slovenia.
  * Xu et al. [2021] Bingke Xu, Yue Cui, Zipeng Sun, Liwei Deng, and Kai Zheng.  Fair Representation Learning in Knowledge-aware Recommendation.  In _2021 IEEE International Conference on Big Knowledge (ICBK)_ , pages 385–392, Auckland, New Zealand, 2021. IEEE.  doi: 10.1109/ICKG52313.2021.00058.
  * Yao and Huang [2017] Sirui Yao and Bert Huang.  Beyond parity: Fairness objectives for collaborative filtering.  In Isabelle Guyon, Ulrike von Luxburg, Samy Bengio, Hanna M. Wallach, Rob Fergus, S. V. N. Vishwanathan, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , pages 2921–2930, Long Beach, CA, USA, 2017. NeurIPs.  URL https://proceedings.neurips.cc/paper/2017/hash/e6384711491713d29bc63fc5eeb5ba4f-Abstract.html.