# Sedimentation of Fractal Aggregates in Shear-Thinning Fluids
## 2 Mathematical Model and Numerical Method
### Numerical Method  
The calculation of the regime drag correction coefficient requires the knowledge of the orbit followed by the aggregate while pulled by the force. This would need, at each time step, the calculation of the particle angular velocity from the solution of the set of equations presented in the previous section (as discussed above, the translational dynamics is irrelevant). Instead of directly solving the governing equations to compute the orbit, we can greatly speed-up the calculations by precomputing a database of particle translational and angular velocities for different orientations of the aggregate [35; 36]. The velocities needed at the right-hand side of Equations (10) and (11) are, then, obtained by interpolating the data of the database. In this way, the computational effort is only due to the construction of the database that can be used to compute the orbit for any initial orientation. Of course, the database must be re-computed for every particle morphology. Due to the irregular particle shape, two orthogonal vectors fixed with the particle are needed to track the evolution of its orientation. However, it should be noted that any configuration obtained by rotating the aggregate around the applied force is equivalent in terms of particle dynamics and drag coefficient. Consequently, one orientation vector is sufficient to fully describe all the possible configurations of the aggregate with respect to the applied force. This can be readily seen if we consider the particle fixed in the laboratory frame and the applied force is rotated. All the possible configurations are, indeed, obtained by rotating the applied force around the unit sphere (so just considering one orientation vector). This is, in fact, the procedure we adopted to build the database, i.e., we fix the aggregate orientation as generated by the particle-cluster algorithm and solve the fluid dynamics problem for several orientations of the applied force uniformly distributed over the unit sphere. More specifically, we divide the unit sphere in a triangular mesh with icosahedral symmetry. The orientations of the force are taken as the directions connecting the center of the unit sphere and the vertices of the icosahedral mesh. We select an icosahedral subdivision with 162 vertices, verifying that this subdivision is sufficient to assure a good accuracy of the interpolation. Indeed, by computing the interpolating functions on an icosahedral grid with 42 vertices, the maximum relative error is \(\sim\)3%. Once the database of particle velocities (and the corresponding drag correction coefficients) has been computed, the rotational dynamics of the aggregate can be calculated by integrating Equation (11). It is, however, more convenient to compute the particle dynamics in the body reference frame, i.e., by fixing the orientation of the particle and rotating the applied force. The adopted procedure is summarized in Algorithm 1. We have denoted with \(f\) the (time-dependent) applied force vector acting on the aggregate, and with \(f_{0}\) its initial value. In step 3, the search procedure described in [37] is used. A linear interpolation over the spherical triangle mesh is done in step 4 [37; 38]. The update of the force \(f\) is carried out through quaternions [39]. Specifically, in step 5, a third-order Adams-Bashforth scheme is used to integrate in time the quaternions through the angular velocity computed at the previous step. A rotation matrix in terms of quaternions is constructed and used to update \(f\). The procedure just described gives the time-evolution of the applied force \(f\) and of the drag correction coefficient \(X(t)\) that can be used to calculate the regime drag correction coefficient \(X_{\mathrm{R}}\). Finally, the aggregate orbit can be reconstructed by transforming the force from the body to the lab reference frame, i.e., by rotating the aggregate according to the rotation matrix that, at every instant, transform \(\mathbf{f}\) in \(\mathbf{F}=(F,0,0)\) (that is the applied force in the lab frame).