... ... @@ -32,7 +32,7 @@ v^{i+1} &= \mathrm{Contacts}\left(x^i, v^i, \vfree^{i+1}, θ\right)\\ x^{i+1} &= x^i + ((1-θ)\,v^i + θ\,v^{i+1})\,Δt \end{align*} où les contacts sont résolus de manière à ce qu'il n'y ait pas d'interpénétration dans $x^{i+1}$. Jusqu'ici, j'ai pris $θ = 0$ mais passer à $θ = 0.5$ ne devrait pas poser de problème. Jusqu'ici, j'ai pris $θ = 1$ mais passer à $θ = 0.5$ ne devrait pas poser de problème. Ce schéma fonctionne mais on est limité à des petits pas de temps (autant pour la précision que pour la stabilité semble-t-il). Afin de s'affranchir de cette limite, un schéma multi-pas partiellement implicite (IMEX RK) sera utilisé pour calculer $\vfree$. ... ...
 \documentclass{paper} \usepackage[a4paper,margin=3cm]{geometry} \usepackage[utf8]{inputenc} \usepackage{amsmath} \def\labelitemi{--} \begin{document} \section{Contacts resolution} \subsection{Potential contacts detection} We start by generating a list of the potential contacts that could occur during a time step. A contact occurs between two objects if they intersect each other before the end of the time step. The potential contacts are all the pairs of objects separated by a distance smaller than an alert distance $A$ at the beginning of the time step. This alert distance should be chosen so that $A > 2\,\Delta t\, v_\text{max}$ where $v_\text{max}$ is the speed of the fasted object. To avoid a large number of potential contacts (which would slow down the resolution), it is recommended to choose $A < r$, r being the radius of the particles, and to compute $\Delta t$ accordingly. A contact $c$ between the objects $o_{c0}$ and $o_{c1}$ is characterized by the distance between those two objects at the beginning of the time step, $D_c$, and a direction $\mathbf n_c$. The relative velocity absorbed by the contact in this direction is noted $dv_c$. The factors $a_{c0}$ and $a_{c1}$ determine how this relative velocity loss is shared by the two objects. The effect of the contact is to remove $a_{c0}\, dv_c\, \mathbf n_c$ from $\mathbf v_{o_{c0}}$ and to add $a_{c1}\, dv_c\, \mathbf n_c$ to $\mathbf v_{o_{c1}}$. To avoid a number of operations in scaling as as the square of the number of particles, we use an octree to limit the search of contacts to the pair of objects for which the distance between their bounding boxes smaller than $A$. \subsubsection* {particle-particle contacts} We only consider spherical particles. The distance between the particles $D_c$ and the contact direction $\mathbf n_c$ are given by: $D_c = ||\mathbf x_{o_{c1}} - \mathbf x_{o_{c0}}|| - (r_{o_{c0}} + r_{o_{c1}}) \text{ and } \mathbf{n}_c = \frac{\mathbf x_{o_{c1}} - \mathbf x_{o_{c0}}}{||\mathbf x_{o_{c1}} - \mathbf x_{o_{c0}}||}.$ where $\mathbf x_i$ is the position of the center of the particle $i$ and $r_i$ its radius. All contacts with $D_c > A$ are discarded. The velocity absorbed by the contact is shared by both particles according to their masses $m_i$: $a_{c0} = \frac{m_{o_{c1}}}{m_{o_{c0}} + m_{o_{c1}}} \text{ and } a_{c1} = \frac{m_{o_{c0}}}{m_{o_{c0}} + m_{o_{c1}}}.$ \subsubsection* {particle-boundary contacts} Three type of boundaries are implemented: sphere, segment and triangle, this is enough to represent any meshed boundary. Our implementation of the contacts between a particle and a segment (resp. triangle) does not handle the intersection with the extremities of the segment (resp. edges of the triangle). For a meshed boundary composed of triangles, it is necessary to consider not only the triangles but also the edges and the vertices (as sphere of radius $0$). We assume that $o_{c0}$ is the boundary and $o_{c1}$ is the particle. The movement of the boundary is always prescribed and is not affected by the contact, so that $a_{c0} = 0$ and $a_{c1} = 1$ for all particle-boundary contacts. The distance $D_c$ and direction $\mathbf n_c$ of the particle-spheres contacts are obtained in the same way as the particle-particle contacts. For the particle-segment (resp. triangle) contacts, we define the point $\mathbf x_p$ as the projection of the center of the particle on the line (resp. plane) containing the segment (resp. triangle), which allows us to define the distance and direction of the contact: $D_c = ||\mathbf x_{o_{c1}} - \mathbf x_p|| - r_{o_{c1}} \text{ and } \mathbf{n}_c = \frac{\mathbf x_{o_{c1}} - \mathbf x_p}{||\mathbf x_{o_{c1}} - \mathbf x_p||}.$ Two cases can lead to a negative $D_c$. If $x_p$ is inside the segment (resp. triangle), there is effectively a (small) intersection at the beginning of the time step, due to round-off or tolerance errors at a previous time step, in this case we set $D_c$ to $0$ and we keep the contact (it is necessary to set $D_c$ to $0$, otherwise we induce a velocity that moves the two objects away from each others at the following time steps by inertia). If $x_p$ is not inside the segment (resp. triangle), the particle trajectory will cross the extremities of the segment (resp. edges of the triangle) before crossing the segment (resp. triangle). The contact with the segment (resp. triangle) is discarded. \subsection{Iterative algorithm} The iterative algorithm to find a solution that satisfies simultaneously all the potential contacts is very simple. We just iterate over all the contacts repetitively until a stop criterion is satisfied. For each contact $c$, we compute $\Delta C_c = C(\mathbf{v}_{\text{locfree}}) - dv^i_c$ where $\mathbf{v}_{\text{locfree}} = (\mathbf{v}^i_{o_{c0}} - \mathbf{v}^i_{n_{c1}}) + dv_c^i\, \mathbf{n}_c$, and we update the velocities accordingly: \begin{align*} dv_c^{i+1} &= dv_c^i + \alpha\, \Delta C_c\\ \mathbf v_{o_{c0}}^{i+1} &= \mathbf v_{o_{c0}}^i + \alpha\, \Delta C_c\, a_{c0}\\ \mathbf v_{o_{c1}}^{i+1} &= \mathbf v_{o_{c1}}^i - \alpha\, \Delta C_c\, a_{c1}. \end{align*} We iterate over $i$ while $r_\text{max} = \max_c(|\Delta C_c|) \leq \varepsilon$. The choice of the relaxation factor $\alpha$ may be problematic. If $\alpha$ is too big, this Jacobi algorithm is unstable, if $\alpha$ is too small, it becomes very slow. In practice, $.33$ in 2d and $.125$ in 3d seems reasonable if all the particles have the same size (i.e. the maximum number of contacts affecting one particle is limited). For particles of various sizes, one could try (I have not done it) to use a different $\alpha$ for each contact, e.g. $\alpha_c = \min(n_0^{-1}, n_1^{-1})$ where $n_0$ and $n_1$ are the number of potential contacts detected for the particles involved in this contact $c$. Another possibility (which is the default behaviour in our implementation) is to use a Gauss-Seidel approach. Inside one single iteration of the iterative solver, we immediately use the new velocities updated by the previous contacts. It has the advantage to be faster and stable for $\alpha = 1$ but it is more difficult to implement in parallel . \subsection{Single contact resolution $C(\mathbf v_\textnormal{locfree})$} $v_n = \mathbf{v}_{\text{locfree}}\cdot \mathbf{n}$ $C(\mathbf{v}_\text{locfree}) = \begin{cases} v_n - \frac{D}{\Delta t}, & \text{if the contact is active}\\ 0, & \text{otherwise} \end{cases}$ A contact is said active'' if, ignoring this contact, an intersection of the two objects would occur before the end of the time step. A particle-particle or particle-sphere contact is active iff $v_n\, \Delta t > D$. For a particle-segment (resp. particle-triangle) contact, the condition $v_n\, \Delta t > D$ only ensures that the sphere intersects the line containing the segment $S$ (resp. the plane containing the triangle $T$) during this time step. The contact with the line (reps. plane) occurs in $\mathbf{x}_n = \mathbf{x}_1 + r_1\,\mathbf{n} + \frac{D}{v_n} \mathbf{v}_{\text{locfree}},$ where $\mathbf{x}_1$ is the center of the sphere and $r_1$ its radius. The contact is active iff $v_n\,\Delta t > D$ and $\mathbf{x}_n \in S$ (resp. $\mathbf{x}_n \in T$). This ignores the contacts with the extremities of the segment (resp. edges and vertices of the triangles). In practice, to cope with the round-off errors, we check that the projection of $\mathbf{x}_n$ on the line $\mathbf{x}_0\mathbf{x}_1$ is inside the segment $]{\mathbf{x}_0\mathbf{x}_1}[$, i.e. Similarly, in the case of a triangle, we check that the projection of $\mathbf x_n$ on the plane $\mathbf x_0 \mathbf x_1 \mathbf x_2$ is inside the triangle $]{\mathbf x_0\mathbf x_1\mathbf x_2}[$. \end{document}