Auto Contacts

Choice of a Formulation for Auto Impact

In the case of a surface impacting on itself, it is possible to use one of the previous formulations if considering certain specificities of the auto-contact. The choice of a formulation will depend on two essential criteria: the quality of the description of the contact and the robustness of the formulation. The selected formulation has to provide results that are as precise as possible in a normal operational situation, while still working in a satisfying way in extreme situations. The node-to-node method provides the best robustness, but the quality of the description is not sufficient enough to simulate in a realistic way those contacts occurring during the buckling of a structure.

The node-to-surface contact is the best compromise. However, it has some limitations, the main one being that it cannot detect contacts occurring on the edges of a segment. The most critical situation occurs when this leads to the locking of a part of one surface to another. This phenomenon, being irreversible, might create irrelevant behavior during the deployment of a structure (example, airbag deployment). An example of a locking situation is shown in the above images. To correct this, it is possible to associate a node-to-surface formulation to an edge-to-edge formulation.

Figure 4.

During the impact of one part on another, it is possible to predict which node will impact on a certain segment. In the case of the buckling of a shell structure, such as a tube, it is impossible to predict where different contacts will occur. It is thus necessary to have a fairly general and powerful algorithm that is able to search for impact candidates.

The detail of the formulation of an algorithm, able to search for impact candidates, will depend on the choice of the contact formulation described in the previous chapter. In a node-to-node contact formulation, it is necessary to find for each node the closest node, whose distance is lower than a certain value. In the case of the edge-to-edge formulation, the search for neighboring entities concerns the edges and not the nodes. However, we should note that in some algorithms, the search for neighboring edges or segments is obtained by a node proximity calculation. Moreover, an algorithm designed to search for proximity of nodes can be adapted in order to transform it into a search for proximity of segments or even for a mixed proximity of nodes and segments.

When using direct search, at each cycle the distance is calculated from each entity (node, segment, edge) to all others. The quadratic cost (N*N) of this algorithm makes it unusable in case of auto-contact.

Topologically Limited Search

In a simulation in fast dynamics, geometrical modifications of the structure are not very important during one cycle of calculation. It is then possible to consider neighborhood search algorithms using the information of the previous cycle of computation. If for a node the nearest segment is known at the previous cycle, it is then possible to limit the search for this node to the segments topologically close to the previous one (the segments having at least one common node). Furthermore, if an algorithm based on the search of neighboring nodes is used, then the search may be limited to the nodes of the segments connected to the previous closest node.

Figure 5.

It remains however necessary to do a complete search before the first cycle of calculation. We will also see that this algorithm presents significant restrictions that limit its use to sliding surfaces (it cannot be used for auto-impacting surfaces).

The cost of this algorithm is linear (N), except at the first cycle of computation, during which it is quadratic (N*N). The combination of this algorithm with one of the following two is also possible.

Bucket Sort

Sorting by boxes consists in dividing space in to steady boxes (not necessarily identical) in which the nodes are placed. The search for closest nodes is limited to one box and the twenty-six neighboring boxes. The cost of this sorting is linear (N) for regular meshes. For irregular meshes, an adaptation is possible but its interest becomes less interesting compared with the next solution. This three-dimensional sorting is of the same kind as one-way sorting with direct addressing or needle sort.

Figure 6.

In order to limit the memory space needed, we first count the number of nodes in each box, thereby making it possible to limit the filling of those boxes that are not empty. In the two-dimensional example shown above, nodes are arranged in the boxes as described in the following table. With this arrangement, the calculation of distances between nodes of the same box is not a problem. On the other hand, taking into account the nodes of neighboring boxes is not straightforward, especially in the horizontal direction (if the arrangement is first made vertically, as in this case). One solution is to consider three columns of boxes at a time. Another solution, more powerful but using more memory, would be for each box to contain the nodes already located there plus those belonging to neighboring boxes. This is shown in the third series of columns in the following table. Once this sorting has been performed, the last step is to calculate the distances between the different nodes of a box, followed by the distance between these nodes and those of the neighboring box. In box 0,3 for example, fifteen distances must be calculated.

11-14, 11-3, 11-8, 11-12, 11-15, 11-16, 11-17, 11-18, 14-3, 14-8,

14-12, 14-15, 14-16, 14-17, 14-18.

In this example you have considered a search based on nodal proximity. It is possible to adapt this algorithm in order to arrange segments or even edges in the boxes. It will be necessary however to reserve more memory, for a segment can overlap several boxes.

Quick Sort

The octree is a three-dimensional adaptation of fast sorting. Space is divided into eight boxes, each one being subdivided into eight boxes. In this way a tree with eight branches per node is obtained. An alternative to the "octree", closer to the quick sort, consists in successively dividing space in two equal parts, according to directions X, Y, or Z. This operation is renewed for each of the two resulting parts as long as some segments or nodes are found in the space concerned. The main advantage to this algorithm, as compared to sorting by boxes, lies in the fact that its performance is affected neither by the irregularity of the mesh, nor by the irregularity of the model. The cost of this algorithm is logarithmic (N*Log(N)).

Figure 7. Possible Impact of Node n₁

In order to illustrate the 3D quick sort, consider a search for node-segment proximity. At each step, space is successively divided into two equal parts, according to directions X, Y or Z. The group of nodes is thus separated into two subsets. We thereby obtain a tree organization with two branches per node. After each division, a check is made to determine whether the first of the two boxes must be divided. If so, it will be divided similar to the previous one. If not, the next branch is then checked. This recursive algorithm is identical to the regular fast sorting one.

The segments are also sorted using the spatial pivot. The result of the test can lead to three possibilities: the segment is on the left side of the pivot, on the right side of the pivot or astride the pivot. In the first two cases, the segments are treated similar to the nodes, but in the third situation, the segment is duplicated and placed on both sides.

The formulation used in Radioss is a penalty type formulation. The choice of the penalty factor is a major aspect of this method. In order to respect kinematic continuity, the penalty spring must be as rigid as possible. If the impedance of the interface becomes higher than those of the structures in contact, some numerical rebounds (high frequency) can occur. To ensure the stability of the integration diagram, without having additional constraints, this rigidity must be low. With a too low penalty, the penetration of the nodes becomes too strong and the geometrical continuity is no longer ensured.

The compromise selected consists in using a stiffness of the same order of magnitude than the stiffness of the elements in contact. This stiffness is nonlinear and increases with the penetration, so that a node is not allowed to cross the surface.

These choices provide a clear representation of physics, without numerical generation of noise, but require the contact stiffness in the calculation of the criteria of stability of the explicit scheme to be taken into account.

After identifying the candidates for the impact, it is necessary to determine whether contact takes place and its precise localization. If for a formulation of node to node contact the detection of the contact is quite easy, it becomes more complex in the case of a node to segment or edge to edge contact. In the case of edge-to-edge contacts, a direct solution is possible if the segments are planar. If not, it is better to triangulate one of the segments, which would then turn it into a node-to-segment contact problem.

The search algorithm for candidates is uncoupled from the rest of the processing of the interfaces. This is not the case with regard to the detection, localization and processing of the contact. These last three tasks significantly overlap with each other so we will limit ourselves to the processing of the contact by penalty for simplicity.

In the case of contact between two solid bodies modeled with 3D finite elements, contact can only take place on the segments of the external surface. This external surface has a certain orientation and the impact of a node can come only from the outside. Most of the node-to-surface contact algorithms use this orientation to simplify detection of contact. In the case of impact of a two-dimensional structure modeled with shell or plate finite elements, contact is possible on both sides of the surface.

For an oriented surface, it is necessary to consider contacts of the positive dimension side of the surface on itself, contacts of the negative dimension side on itself, and the contact of the "positive" side on the "negative" one. This last situation, which is quite rare, can occur in the case of a surface rolled up into itself or during the impact of a twisted surface.

Node to Segment Contact

The use of a "gap" surrounding the segment is one way of providing physical thickness to the surface and makes it possible to distinguish the impacts on the top or on the lower part of the segment. The contact is activated if the node penetrates within the gap or if the distance from the node to the segment becomes smaller than the gap.

To calculate the distance from the node to the segment, we make a projection of this node on the segment and measure the distance between the node and the projected point.

Figure 8.

The projected point is calculated using isoparametric coordinates for a quadrangular segment and barycentric coordinates for a triangular segment. In the case of any quadrangular segment, the exact calculation of these coordinates leads to a system of two quadratic equations that can be solved in an iterative way. The division of the quadrangular segment into four triangular segments makes it possible to work with a barycentric coordinate system and gives equations that can be solved in an analytical way.

From the isoparametric coordinates ($$s$$, $$t$$) of the projected point ($$Ni$$), we have all the necessary information for the detection and the processing of contact. The relations needed for the determination of $$s$$ and $$t$$ are: the vector $$NiNs$$ is normal to the segment at the point $$Ni$$; and the normal to the segment is given by the vectorial product of the vectors $$si$$ and $$ti$$.

$\begin{array}{l}H1=(1-s)(1-t)/4\\ H2=(1+s)(1-t)/4\\ H3=(1+s)(1+t)/4\\ H4=(1-s)(1+t)/4\\ \overrightarrow{ON}i=H1\overrightarrow{ON}1+H2\overrightarrow{ON}2+H3\overrightarrow{ON}3+H4\overrightarrow{ON}4\\ \overrightarrow{s}i=(1-t)\overrightarrow{N1N2}+(1+t)\overrightarrow{N4N3}\\ \overrightarrow{t}i=(1-s)\overrightarrow{N1N4}+(1+s)\overrightarrow{N2N3}\\ \overrightarrow{n}i=(\overrightarrow{s}i\times \overrightarrow{t}i)/(\overline{\overrightarrow{s}i\times \overrightarrow{t}i})\\ \overrightarrow{NiNs}=a\overrightarrow{n}i\end{array}$

After bounding the isoparametric coordinates between +1 or -1, the distance from the node to the segment and the penetration are calculated:

$\begin{array}{l}D=\overline{NiNs}\\ P=\max (0,Gap-D)\end{array}$

A penalty force is deducted from this. It is applied to the node $$Ns$$ and distributed between the four nodes ($$N1$$, $$N2$$, $$N3$$, $$N4$$) of the segment according to the following shape functions.

Edge to Edge Contact

The formulation of edge-to-edge contact is similar to that of node-to-segment contact. The gap surrounding each edge defines a cylindrical volume. The contact is detected if the distance between the two edges is smaller than the sum of the gaps of the two edges. The distance is then calculated as:

$\begin{array}{l}2\overrightarrow{N1Ni}=(1-s)\overrightarrow{N1N3}+(1+s)\overrightarrow{N1N4}\\ D=\left|\overrightarrow{N1Ni}-\left(\overrightarrow{N1N}i\overrightarrow{N1N2}/{\overline{N1N2}}^2\right)\overrightarrow{N1N2}\right|\\ \partial D/\partial s=0\end{array}$

The force of penalty is calculated as in node-to-segment contact. It is applied to the nodes $$N1$$, $$N2$$, $$N3$$, $$N4$$ and therefore ensures the equilibrium of forces and moments.

Figure 9.