In this section, numerical procedures and implementation-specific aspects related to particle tracing are discussed.
Given a velocity field exactly defined at any time at any point in the model space, and given the injection (emitter) point and time of a particle, one likes to integrate the velocity field to contain the trajectory (or trace) of that particle. Since the instant, at which a point along the trajectory is passed by the particle, is defined, one can can make an exact path integration of the particle. However, in a CFD calculation, the velocity field is not exactly defined but interpolated in space and time, thus the accuracy of the integration deteriorates with growing path length. This has to be taken into account when interpreting simulated particle traces. This section discusses how baspl++ performs this integration.
In CFD applications, the velocity fields are calculated at discrete instants. To obtain the velocity field at any instant, the values of the calculate fields are linearly piece-wise interpolated in baspl++.
With these values, a fourth-order integration (Runge-Kutta) of the time-dependent velocity field is performed to obtain the next point of the particle's trace. This integration scheme has shown to be sufficiently precise, as other sources of numerical errors are more dominant. As with all numerical integration methods, the choice of the step size is very important. Many applications perform an adaptive algorithm, where the step size is varied until a good performance/error ratio is achieved. We found that step size adaptation is costly as it always requires at least two times more operations than would be required if the optimal step size was known in advance. Therefore and for the reasons outlined in the next paragraph, we decided to use for baspl++ a constant step size, constant with respect to each element, so that passing through an element requires approximately 5 integration steps.
This number of steps is sufficient for good accuracy, as the velocity field inside a linear element (and baspl++ handles linear elements only for particle tracing) cannot change rapidly. For hexahedral elements, even if the velocity field at the vertices point in very different directions, the angle between consecutive trajectory points remains below 15 degrees. In unstructured meshes, where tetrahedral elements predominate, the particle trace often passes near a vertex of the tetrahedron. In this case it makes no sense either to choose a smaller step size, as the value of the velocity field (due to the geometrically linear interpolation) is almost entirely defined by the value at that vertex.
To be able to set the step size in terms of the element's size, one is forced to integrate over the velocity field in element-internal coordinates. That is, the global velocity field, which is defined at the element's vertices, is transformed into internal coordinates, so that the velocity vectors at the vertices are in internal coordinates. This transformation is achieved by solving a linear system of the Jacobian matrix at the current position, the internal velocity and the global velocity. Integrating then yields the new position in element-internal coordinates. An important benefit of this approach is that calculating the position in internal coordinates from global coordinates is only necessary when passing from one element to another.
Of interest is also the velocity interpolation inside elements, as its quality influences the accuracy of the integration process. Many post-processors split up hexahedral, prism, and pyramidal elements into tetrahedrons. The reasoning behind is a gain in computational efficiency, since the Jacobian matrix of tetrahedral elements is constant over the whole element volume and thus has to be calculated only once. However, with such an approach the geometrically non-linear shape functions of e.g. hexahedral elements are replaced by piece-wise linear functions, thus the accuracy of the interpolation is impaired, especially in turbulent regions. We decided to use the respective Finite Element shape functions for first-order hexahedrons, pyramids, prims/wedges, and tetrahedrons to interpolate the velocity field. The possible loss of speed is outweighed by the higher fidelity of the calculated traces. baspl++ goes to great lengths to calculate the Jacobian for the different element types (tetrahedron, pyramid, prism, hexhadron) in a most efficient manner, thus the gap between a tetrahedral approach and baspl++'s is smaller than one would normally expect.
There are structured finite volume meshes where the values are defined per-cell and not on the nodes of the mesh. baspl++ constructs a new mesh so that nodes coincide with cell centers. To make the new mesh structured too, missing nodes on the mesh boundary are interpolated. The mesh is then treated like any ordinary hexahedral mesh, with values defined on the nodes, thus trilinear interpolation inside hexahedral cells is used. An inconvenient of this approach is that at the boundaries of different branches, where - for the interpolation of the missing nodes - only the given mesh is considered and not neighbouring meshes. In many cases this leads to small holes where at least three subdomains join. baspl++'s tracing mechanism has special provisions for the case where a particle trace falls into such a hole, as is described in the next paragraph.
Finally, the passing from one element to the next can pose problems as well. baspl++ maintains a table where the element-element neighbourhood relationships are recorded (neighbouring elements have an element face in common). This table is built once when the tracing facility is started and provides a very fast element lookup. But sometimes, depending on the mesh type, the exact neighbour does not exist (e.g. Chimera meshes). A special element octree serves in this case as a fallback mechanism to find the element. In case there are "holes" in the mesh, that is, the particle is not surrounded by an element, the velocity field is kept constant and the particle moves until it is again inside an element. Thus baspl++ can handle any type of mesh type or combination of different mesh types.
