The electrostatics were introduced to the simulation by allowing particles to experience triboelectric charging. The triboelectric charging was based on an effective work function model that characterized triboelectric charging tendency of each material by an effective work function \(\phi\).
The charge transfer rate during a collision is given by
\[\frac{dq}{dt} = H \left(\frac{dA}{dt}\right) \frac{dA}{dt} \frac{\varepsilon}{\delta_c e}\left( \Delta \phi + \mathbf{E}_{ij} \cdot \mathbf{n}_{ij} \delta_c e \right), \label{eq:t1}\]
where \(A\) is the contact area, \(H(\cdot)\) is Heaviside function, \(\delta_c\) is the maximum distance where the charge can transfer from object to object, typically taken to be \(500 \nano \meter\), \(e\) is the charge of an electron, \(\varepsilon\) is the electric permittivity, and \(\mathbf{E}_{ij}\) is the electric field at the contact point.
The Heaviside function is used to stop the charge transfer back to where it came from. It is generally accepted that there is some level of charge relaxation happening when particles are coming apart that is neglected in the Eq. \eqref{eq:t1}. However, there is no prior work or model to determine this backflow of charge, and it will be neglected in this study. Therefore, the Eq. \eqref{eq:t1} will slightly overestimate the charge transfer between the particles.
The electron field at the contact point was computed from the neighboring particles that included also the colliding particles by direct sum
\[\mathbf{E}_{ij} = -\frac{1}{4 \pi \varepsilon} \sum_{k} \frac{ q_k \mathbf{r}_k }{ \| \mathbf{r}_k \|^3}, \label{eq:t2}\]
where \(\mathbf{r}_k\) is vector pointing from \(k\):th particles center to the center of the contact. The neighboring particles are defined by setting a cut-off radius for the electrostatic interactions. Any particles that are further than the predefined cut-off radius are neglected from the sum.
During the simulation the charge transfered between the particles \(i\) and \(j\) between timesteps \(n\) and \(n+1\) becomes
\[\Delta q^{n+1} = H( \Delta A ) \Delta A \frac{\varepsilon}{\delta_c e}\left( \Delta \phi - \mathbf{E}^n_{ij} \cdot \mathbf{n}_{ij} \delta_c e \right),\]
where the electric field \(\mathbf{E}^n_{ij}\) is evaluated at the previous time step \(n\), and \(\Delta A\) is the change of contact area between the timesteps. The transfered charge \(\Delta q^{n+1}\) is typically small, and the stiffness of the Eq. \eqref{eq:t1} does not become an issue with timesteps typically used for the soft sphere model. The good sides of the triboelectric model of Eq. \eqref{eq:t1} is that it can capture collisions of different materials, the effect of surrounding electric field, and the effect of contact force via the contact area. The subtle physics can be introduced by making the effective work function \(\phi\) as a function of various parameters such as particle size or ambient humidity.
The electrostatic force on a particle was evaluated as
\[\mathbf{F}_i = q_i \mathbf{E}(\mathbf{x}_i), \label{eq:t3}\]
where \(q_i\) is the charge of the particle, \(\mathbf{x}_i\) is the location of the \(i\):th particle, and \(\mathbf{E}(\mathbf{x}_i)\) is the electric field at the particles location. The electric field is evaluated from neighboring particles as before by the pairwise summation
\[\mathbf{E}(\mathbf{x}_i) = -\frac{1}{4 \pi \varepsilon} \sum_{k \neq i} \frac{q_k (\mathbf{x}_k - \mathbf{x}_i )}{ \|\mathbf{x}_k - \mathbf{x}_i\|^3}. \label{eq:t4}\]
In our simulations each particle presented individual particles and not collection of particles like in parcel simulations, hence no screening is needed in Eq. \eqref{eq:t4}.
Recall that the charge transfer of a particle during a collision is given by
\[\frac{dq}{dt} = H\left(\frac{dA}{dt}\right) \frac{dA}{dt} \frac{\varepsilon}{\delta_c e}\left( \Delta \phi - \mathbf{E}_{ij} \cdot \mathbf{n}_{ij} \delta_c e \right), \label{eq:e1}\]
where \(\mathbf{E}_{ij}\) is electric field in the contact point, and \(\Delta \phi\) is the workfunction difference. The electric field in the contact point will depend on the particles charge. Therefore, the right hand side of Eq. \eqref{eq:e1} will depend on the particles charge. The Eq. \eqref{eq:e1} can be written in terms of particles charge transfered by splitting the charge to initial charge \(q_0\) and denoting the transfered charge by \(q\). In particular, we have \(q(0)=0\) as at the beginning of contact no charge is transfered. This allows us to write
\[\frac{dq(t)}{dt} = \frac{dA}{dt} \left( \beta - \alpha q(t) \right), \label{eq:e2}\]
where \(\alpha\) is a qeometrical coefficient depending only on particle radius, and \(\beta\) is coefficient depending on the collision type (wall or other particle) and initial charge(s). The detailed formulas for \(\alpha\) and \(\beta\) are given in the appendix. Solving the ODE \eqref{eq:e2} for \(q\) with initial condition \(q(0)=0\) and \(A(0) = 0\), and assuming that the contact area is independent function of charge \(q(t)\) yields
\[q(t) = \alpha\beta \left( 1 - e^{-\alpha A(t)} \right). \label{eq:e3}\]
Let \(t=T\) be time where contact area \(A(t)\) obtains its maximum value \(A_{max}\). Inserting now \(t=T\) to Eq. \eqref{eq:e3} yields charge transfered during one collision
\[\Delta q = \alpha\beta \left( 1 - e^{-\alpha A_{max}} \right). \label{eq:e4}\]
For typical collisions in a fluidized bed we have \(A_{max} \ll \pi r_{eff}^2\) where \(r_{eff}\) is the effective radius of the Hertzian model. This yields \(\| -\alpha A_{max} \| \ll 1\). Proof of this is given in the appendix. Invoking Taylors series expansion for \( \exp(-\alpha A_{max})\) and neglecting higher order terms gives
\[\Delta q = \alpha^2 \beta A_{max}. \label{eq:e5}\]
The ratio between charge transfer of a real and simulated soft particle is therefore,
\[\frac{ \Delta q^{real} }{\Delta q^{soft}} = \frac{ \beta^{real} A_{max}^{real} }{ \beta^{soft} A_{max}^{soft} }.\]
To tackle the problem of soft particles againts the real hard particles we assume for time being that the maximum area is independent on the coefficient of restitution with. Setting coefficient of restitution to one allows us to inspect the energy balances of the collision. The potential energy of Hertzian spring is given by
\[E_{p} = \int F d\delta = \frac{8}{15} Y_{eff} \sqrt{ r_{eff} } \delta^{5/2}, \label{eq:e6}\]
where \(Y_{eff}\) is the effective Youngs modulus, and \(\delta\) is the overlap distance. The overlap distance is connected to contact area by
\[A = \pi r_{eff} \delta.\]
When the contact area obtains its maximum value we have
\[\frac{dA}{dt} = \pi r_{eff} \frac{d\delta}{dt} = \pi r_{eff} v_{\bot} = 0,\]
where \(v_{\bot}\) is the normal velocity. Hence at the maximum contact area we have \(v_{\bot} = 0\), and the normal kinetic energy is transformed totally to Hertzian springs potential energy and to potential energy of the electric field. Assuming that the potential energy of the electricfield is negligible to the springs potential at the maximum contact area energy yields \(E_p = \frac{1}{2} m v_{\bot}^2\).
Solving for \(\delta_{max}\) from Eq. \eqref{eq:e6} we obtain
\[\frac{A_{max}^{real} }{ A_{max}^{soft} } = \left( \frac{Y^{soft}}{Y^{real}} \right)^{2/5} \left( \frac{v_{\bot}^{real} }{ v_{\bot}^{simu} } \right)^2, \label{eq:e7}\]
where \(Y\) stands for particles Youngs modulus. The initial collision velocities may change slightly due to Youngs modulus, but are assumed to depend mainly on the hydrodynamics and the overal charge, hence we have \(v_{\bot}^{real} \approx v_{\bot}^{simu}\). Hence, to correct the charge transfer in a single collision in a soft simulation one has to multiply the charge transfer by
\[\sigma = \frac{ \Delta q^{real} }{\Delta q^{soft}} = \left( \frac{Y^{soft}}{Y^{real}} \right)^{2/5}. \label{eq:ee7}\]
Shifting our attention to particles charge over the total simulation time, and denoting that by \(Q(t)\) we have
\[\frac{dQ(t)}{dt} = \omega(t) \Delta q(t), \label{eq:e8}\]
where \(\omega(t)\) is number of collisions per time unit, and \(\Delta q\) charge transfered during one collision that is obtained from Eq. \eqref{eq:e5}. We assume here that the collision rate would be only function of time and particle parameters even it may depend on the charge in a complicated way.
As noted before \(\Delta Q(t)\) depends on \(\beta(t)\) and \(v_{\bot}\) (now function of time as it depends on \(Q(t)\)) that is an affine function of the particles charge before each collision. Therefore we have
\[\frac{dQ(t)}{dt} = \omega(t) a \left( Q_{eq}-Q(t) \right), \label{eq:e9}\]
where \(a\) is a coefficient that depends on \(v_{\bot}\) and Younds modulus, and \(Q_{eq}\) is coefficient that depend only on the particles initial charge at the beginning of the whole simulation and on the imposed work functions. The exact values of the coefficients are given in the appendix and depend on the collision type. In particular coefficient \(a\) depends linearly on the maximum contact area, and \(Q_{eq}\) is independent on the Youngs modulus.
Solving again the familiar ODE \eqref{eq:e9} yields
\[Q(t) = Q_{eq} \left( 1 - e^{ -a \int \omega(t) dt } \right). \label{eq:e10}\]
Assume that the maximum contact area of soft simulation would be multiplied by a number, say \(\sigma\). Setting the right hand sides of Eq. \eqref{eq:e10} as equal for soft and hard particles, and after some algebraic manipulation we obtain
\[\sigma(t) = \left( \frac{Y^{soft}}{Y^{real}} \right)^{2/5} \left( \frac{v_{\bot}^{real}(t) }{ v_{\bot}^{simu}(t) } \right)^2 \frac{\int_0^{t} \omega_{real}( \tau ) d\tau }{\int_0^{t} \omega_{soft}(\tau) d\tau } \label{eq:e11}\]
The Eq. \eqref{eq:e10} is valid if the charging does not affect the maximum contact area during one collision. The Eq. \eqref{eq:e10} has few implications. Firstly, to simulate a system of hard particles with soft particles it is enough to multiply soft particles contact area by coefficient given in Eq. \eqref{eq:ee7} to obtain correct charging rate. This correction is valid as long as charges in the system are small enough not to alter the collision rate function \(\omega(t)\). If the collision rate is changed the coefficient is given by Eq. \eqref{eq:e11}.
The second implications is that if we have simulated a charging rate of a system with soft particles the results should be shown with time scaled by the coefficient of Eq. \eqref{eq:e11} or by a coefficient of Eq. \eqref{eq:ee7} if the collision dynamics is not changed by the charging.
The proposed correction of Eq. \eqref{eq:ee7} was tested in number of two particle collision with varying coefficients of restitution, Youngs modulus, and various initial velocities. The charge on the particles was fixed as \(\pm 1.0^{-13}\mathrm{C}\), and the work function difference was zero.
[][c][0.75]\(Y=10^5\) [][c][0.75]\(Y=10^6\) [][c][0.75] \(Y=10^7\) [][c][0.75] \(Y=10^8\) [][c][0.75]\(v\) [m/s] [][c][0.75]\(\Delta q\) [C] [][c][0.75]\(e=0.9\)
\label{fig:e1}
[][c][0.75]\(Y=10^5\) [][c][0.75]\(Y=10^6\) [][c][0.75] \(Y=10^7\) [][c][0.75] \(Y=10^8\) [][c][0.75]\(v\) [m/s] [][c][0.75]\(\Delta q\) [C] [][c][0.75]\(e=0.9\)
\label{fig:e2}
[][c][0.75]\(Y=10^5\) [][c][0.75]\(Y=10^6\) [][c][0.75] \(Y=10^7\) [][c][0.75] \(Y=10^8\) [][c][0.75]\(v\) [m/s] [][c][0.75]\(\Delta q\) [C] [][c][0.75]\(e=0.3\)
\label{fig:e3}
The uncorrected case is shown in Fig. \ref{fig:e1} showing very different values of charge transfer. The same case was run again with corrected charge transfer, and the results are shown in Fig. \ref{fig:e2}. The corrected case shows very neat collapse of all the cases to the base case with Youngs modulus \(10^8\). The cases in Fig. \ref{fig:e1} and \ref{fig:e2} are computed with high coefficient of restitution \(0.9\) that is close to coefficient of restitution \(1.0\) assumed in our theoretical reasoning. We also tested the correction in lower coefficient restitution \(0.3\) case to see if the results still showed good agreement. The results are shown in Fig. \ref{fig:e3}, and the curves collapsed fully. In the very small velocity limit, the coefficient of restitution caused the particles to come at full stop. Suprisingly, the corrected charge transfer was still to correct value even the coefficient of restitution would dominate the collision in such case. In addition we also tested multiple particle simulations with wall with different Youngs modulus, namely \(10^5\) and \(5 10^6\), with the correction. The correction was taken respect to Youngs modulus \(3\times 10^9\) that corresponds to polyethylene. The results showed very good fit, and the curves resembled each other exactly. This is not that suprising, since once the charging of each individual collision is captured accurately, the overall physics should remain the same.