ACTA MECHANICA SINICA (English Series), Vo1.16, No.4, November 2000 The Chinese Society of Theoretical and Applied Mechanics Chinese Journal of Mechanics Press, Beijing, China Allerton Press, INC., New York, U.S.A.

A CONTACT SEARCHING ALGORITHM CONTACT-IMPACT PROBLEMS* Wang Fujun ( : t ! ~ $ )

ISSN 0567-7718

FOR

Cheng Jiangang ( ~ t ~ )

Yao Zhenhan (Ngt~N) (Department of Engineering Mechanics, Tsinghua University, Beijin 9 100084, China) A B S T R A C T : A new contact searching algorithm for contact-impact systems is proposed in this paper. In terms of the cell structure and the linked-list, this algorithm solves the problem of sorting and searching contacts in three dimensions by transforming it to a retrieving process from two one-dimensional arrays, and binary searching is no longer required. Using this algorithm, the cost of contact searching is reduced to the order of O(N) instead of O(Nlog2N ) for traditional ones, where N is the node number in the system. Moreover, this algorithm can handle contact systems with arbitrary mesh layouts. Due to the simplicity of this algorithm it can be easily implemented in a dynamic explicit finite element program. Our numerical experimental result shows that this algorithm is reliable and efficient for contact searching of three dimensional systems. K E Y W O R D S ; contact impact, contact search, finite element method

1 INTRODUCTION Investigations of the effects of structural impacts are of great importance in many engineering fields. Considering the potential large deformations, a contact can occur anywhere within the structure. The ability to determine contact areas of a general structural system is a fundamental requirement of the contact-impact analysis [1]. Usually, the contact searching includes two different processes [1~3]. The first is the global search that roughly finds all possible candidate contact nodes of surface segments, the second is the local search that locates exactly contact positions after the global search. This paper is concentrated on the global search. So far, a number of various algorithms have been used for the global search. Of them, a typical one is the bucket-sorting algorithm [4] that is embedded in the master-slave contact algorithm[q and the single-surface contact algorithm [2] of DYNA3D by Hallquist. It performed sorting and searching in three dimensions in a nested manner. Later Oldenburg and Nilsson reduced the three-dimensional search to sort and search only in one dimension, that is the position code algorithm [I]. Probing appropriate position codes in a position code vector is done by the binary searching. Thus the complexity of the search is reduced substantially so that the cost of operations of the global search is O(Nlog2N ). In paper [6], Received 29 January 2000, revised 21 July 2000 * The project supported by the National Natural Science Foundation of China (59875045), and the State Key Laboratory of Automobile Safety and Energy Saving (K9705)

Vol.16, No.4

Wang Fujun et al.: A Contact Searching Algorithm

375

an NBS contact searching algorithm for large-scale discrete element simulations is proposed, with total search cost proportional to N. However, the NBS algorithm can be only applied to systems comprising spherical bodies of similar size [6]. In this paper, we propose a new algorithm for global search. It is derived by using the linked-list [7] that is widely used as a data structure in computer programming, and transforms the problem of sorting and searching in three dimensions to a process of sorting and searching within a one-dimension array. It does not need the binary searching, and reduce the searching cost to the order of O(N). Moreover, this algorithm can handle contact systems without limitation of their mesh layouts. 2

FINITE ELEMENT PROBLEMS

FORMULATIONS

FOR CONTACT-IMPACT

As described by Bathe [s] , the equation of the system motion can be transformed into the following discrete equilibrium equation using the standard finite element procedure

M t ~ =tQ _t F +tEe

(1)

In Eq.(1), M is the mass matrix, t/i is the acceleration vector, tQ is the external load vector, tF is the internal load vector and tFc is the nodal contact force vector, tFc can be calculated according to the result of contact searching using the penalty method [2], the Lagrange multiplier method [9,1~ or other methods. Suppose that M is a diagonal mass matrix, Eq.(1) can efficiently be integrated with central differences. The displacements at time t + At are obtained from their values at time t - At and t t+Atu M -1 [At2(tQ -t F +tEe) + 2Mtu - Mt-Atu] (2) -~-

where At is the time increment. In Eq.(2), the contact force tFc, at time t, is the only unknown variables. To calculate the contact force tF~, the contacting boundaries have to be determined using a contact searching algorithm. 3 NEW CONTACT

SEARCHING

ALGORITHM

The goal of the global search procedure is to find out those nodes that are candidates for contact from all pre-defined contact nodes in the structure. We propose a new algorithm through which the contact searching is accomplished by using a cell structure. The cell structure is fixed in space, and the cell sizes are considered to be close to or larger than the average element sizes. The basis of this searching procedure is to identify all nodes in a cell. The key part of this algorithm is the utilization of a special linked-list. 3.1

Cell S t r u c t u r e a n d L o c a t i o n o f N o d e s To identify these nodes and elements that may contact, the concept of cell [11] is introduced. A typical cell structure is shown in Fig.1 that the cells are of equal size in all directions and their edges are parallel to the x, y and z axis. The cell domain is described as follows (Xmin, Xmax, Nx)

(Ymin,Ymax,-/u

(Zmin,Zmax,Nz)

(3)

376

ACTA MECHANICA SINICA (English Series)

Here, Xmin and Zmax are the minimal and maximal x coordinates in the cell domain and l\~ is the number of cells in the x direction. Other variables are analogous to Xmi., Xm~x and N~. The total cell number is calculated by Eq.(4). Nc~n = N~ x N v x N.~

cells "..

z

2000

_1-..

surfacesegments

N~"

x

(4)

The cell location of a node with coordinates (x, y, z) is computed in two steps. First the integer cell numbers in the x , y and z directions are computed by Eq.(5).

Nyl y

Fig.1 Cells enclosing a surface

I~ -- X ~ ( x - 3Cmin)/(Xmax -- ,Train) / - Yrain)/(Ymax

- Ymin)

fy

Ny(y

I:

N z ( z - Zmi,,)/(Zmax -- Zmi.)

(5)

These numbers are then used to compute the cell number of the node, Ic, as shown in Eq.(6). I c = Iz x N~ • Nv + Iy x Nx + I= (6) The size of cells is presumably close to or larger than the average size of elements, as a comparison with the work done in paper [11] in which the ceils are substantially larger than elements and may include too many elements and nodes.

3.2 Constructing the Linked-List Mapping from a set of nodes, {1, 2, 3, 9.., N}, to a set of cells, {1, 2, 3 , . . . , Nc~u}, is defined in such a way that each node is assigned to a cell by its cell number. Figure 2 presents an example where the number of nodes is 9 and the number of cells is 6. Node 1 is assigned to cell 2, node 2 is to cell 3, node 3 to cell 6, node 4 to cell 3, etc. Nodes are mapped by looping over all

cell6 cell 5 cell4 cell 3 cell 2 cell 1

nodes in the ascending numerical order, Fig.2 Mapping of nodes to cells checking their cell number Ic and pushing them along the lists as shown schematically in Fig.3. Thus, the linked-list is formed. For example, the list L3 for cell 3 in Fig.3 was assembled by placing the node 2 onto the list, and then pushing it by node 4, which is pushed by node 5, which is pushed by node 6, which is pushed by node 9, which is the last node added to the list. The linked-list in Fig.4 is represented by two integer arrays, head and n e x t . The first array head contains the last node number mapped to each cell. The array head is a 1-D array with a size of Nceu, where Ncell is the total number of cells. The second array n e x t is a 1-D array with a size of N, where N is the total number of nodes. For each node, the array n e x t contains the next node number in the singly connected list.

Wang Fujun et al.: A Contact Searching Algorithm

Vol.16, No.4

node number

cell number

nodenumber 49-"----'7~

:E 4I : 3

cell number

9

6

5

377

4

6 5 4 3 2 1

2

2

1

:' 7 ,:' 6

8 -1 -1 1 -1

1

head

Fig.3 Lists of nodes for each ceil

3

'N I

-1 next

Fig.4 Numerical representation of the singly connected list

For both arrays, a negative number is used at the end of a singly connected list. So if there are no nodes in a particular cell, a negative number is assigned to the corresponding element of the array head. In short, h e a d represents a pointer to a singly connected list of nodes for each cell, n e x t points to the next element in the list, and ends with a negative number. For example, cell 1 in Fig.2 has no node mapped to it, so the singly connected list of nodes mapped to this cell list L1 is empty, as shown in Fig.3, which is achieved by setting head[l] = - 1 (in Fig.4). In a similar way, the last node mapped to the cell 3 (the list L3) is node 9, thus head[3] = 9. The next node on this list is 6, thus next[9] = 6; the next node is 5, thus next[6] = 5; the next node is 4, thus next[5] = 4; the last node on this list is 2, thus next[4] = 2 and next[2] = - 1 (see the highlighted terms in Fig.4). 3.3 T h e G l o b a l S e a r c h i n g The global searching examines each segment for the potential contacting nodes. With the linked-list, the contact searching can be done efficiently. All the nodes in the cells intersected by a particular segment are thought to be the candidate contacting nodes. The details of the contact searching procedure for a given segment are as follows. First a segment territory must be constructed. For this purpose, the segment boundaries are defined by E={(z~,x2,

x z J~lxmin < I i - xi < - x max

i = 1,2,3}

(7)

where Xl, x2 and x3 are Cartesian components of a position vector x and xmin

9 1 2 = mln(xi, x i , 9" ' , x N'r i )

Nseg is the number of nodes of the segment, segment. The segment territory is defined as T=

{(Xl,X2,

X \ ~ min

z)lXi

-c<_xi

max

xi

= max(x i1, x i2, - - . , x iNsegX )

(8)

x iI is the ith coordinate of node I of the

<_ x m a X + c ,

i=1,2,3}

(9)

where c is a parameter defining the territory expansion. The global searching does not need to be performed at every time intervals. The number of step intervals between two global search cycles can be evaluated by

ACTA MECHANICA SINICA (English Series)

378

lls

2000

\/)max ~ tmax j

where Vmax is the maximum relative velocity of nodes with respect t;o t.he contact surface, and Atm~x is the maximal time interval. If tile territory of a specified segment has been constructed , the number of cells intersected within tile segment territory can be determined by equations (9), (5) and (6). All nodes which are located in those cells can be retrieved from the linked-list (i.e. the arrays head and next) formed in the previous section. These nodes are matched with the segment as candidate contacting nodes, then the local search process is performed to determine the exact contact points on the surface and the distance (gap or penetration) between every node and the segment.

3.4 Local C o n t a c t Searching Whenever a node is found approaching a surface segment, the more exact checks are made in the local search procedure where the actual testing for contact takes place. There are a number of searching algorithms used for the local contact search [s'"'l~ In this project, the inside-outside algorithm[ al is employed. Detailed discussion on the algorithm is stated in Reference [3]. 3.5

The Implementation of the New Algorithm The contact searching procedure in accordance with the new contact searching algorithm can be summarized in the algorithmic form as follows: (1) loop over all nodes in the contact system

{ 9 calculate the cell number of the current node according to its coordinates 9 place tile node number to the array head and next according to cell number

} (2) loop over all element segments in the contact system

{ 9 determine the segment territory of the current segment 9 find the cells intersected by the segment territory (3) loop over all cells intersected by the segment territory

{ 9 get all nodes in the current cell from array head and next (4) loop over all nodes in the cell

{ 9 check for contact between the current node and the current element segment by local contact search process

} } } 3.6 General Remarks There are 4 loops in the algorithm presented in the section 3.5, which consist of 1 loop over all nodes loop (1), 1 loop over all element segments loop (2), 1 l o o p over cells intersected by one segment t e r r i t o r y - - - - l o o p (3), and 1 loop over nodes in one cell loop (4). Assuming that the number of element segments approximately corresponds

Voi.16, No.4

Wang Fujun et al.: A Contact Searching Algorithm

379

to the number of nodes, N, in the contact system, the loop (1) and loop (2) approximately loop over N, separately. The loop (3) and loop (4) depend on neither the number of cells nor the number of nodes in the contact system. Operations performed inside every loop are independent of either the number of cells, or the number of nodes. Therefore, the CPU time to search all contact is proportional to the total number of nodes, N, i.e. T = k. N

(11)

where T is the total search time, N is the total number of nodes, and k is a constant independent of either number of cells or number of nodes in the contact system. Hence, the cost of contact searching is of the order of O(N). RAM requirement related to contact search is easily calculated as the total space occupied by arrays head[Ncen] and next[N]. Thus M = Ncen + N (integers)

(12)

In this contact searching algorithm, the efficiency and memory requirement are influenced by the choice of two user-supplied parameters, the expansion of segment territory with respect to the segment confines and the size of the cells. With a larger territory expansion, the global searching frequency decreases. However, the number of candidate contact nodes involved in the local search procedure increases. With larger cells, the number of contact nodes that are involved in the test with respect to the segment territory increases, but the memory requirement decreases. This implies that an optimal combination of the territory expansion and the cell size should exist for a specific problem. In general, it is acceptable that the size of cells is close to or larger than the average size of elements and the global search is necessary only every 10 to 20 time steps, by our experiences. For a large-scale system, the memory requirement is relatively large if using cells that has a small size. But this is not a difficulty in a parallel computer environment. 4 NUMERICAL

EXAMPLES

The new contact searching algorithm has been implemented into our dynamic explicit finite element code. The inside-outside local contact searching algorithm[ a] is employed in this study to conduct the local searching. For the satisfaction of the contact c o n d i t i o n s - no interpenetrating allowed a standard penalty approach [2] is used. In this finite element program, four-node degenerated shell elements based on Mindlin-Reissner shell theory are employed. In order to demonstrate the general behavior of the contact searching algorithm, three examples are presented. 4.1

I m p a c t of a S q u a r e - S e c t i o n B e a m onto a Rigid Wall To check the performance of the contact searching algorithm for a single-surface contact situation, impact of a square-section beam onto a rigid wall is considered. The same problem has been analyzed with different contact algorithms by many researchers, such as Hallquist[ 2], Zhong [~ Belytschko [1~ and Bathe [12]. The geometry is shown in Fig.5, where length L is 0.30m, width/height a is 0.06m, thickness t is 0.0015m. The material parameters can he found in Ref.[9]. The beam has an initial velocity of 15.6m/s and is free from any external loading. A lumped mass of 1440 kg is attached to the left end of the beam to simulate crash loading on the beam by

380

ACTA MECHANICA SINICA (English Series)

other components connected to the beam. Due to symmetry only one quarter of the beam is modeled. The finite element model consists of 829 nodes and 756 shell elements. In order to perform the contact search, a 60 x 8 x 8 cell structure is constructed. The size of a cell is close to tile size of an element. The analysis has been performed on a Pentium III 450 computer. The analyzed response time is 12ms. Deformed configurations of the beam are shown in Fig.6. Tile total CPU-time required is 35032s, 14.8% of which, i.e. 5185s, is used by the contact searching algorithm.

2000 rigid wall

Attached Mass, M

beam ~V L section

Fig.5 Impact of a square-section beam onto a rigid wall

(a) t = 6 m s (b) t = 1 2 m s Fig.6 The deformed configurations 4.2

Impact Between Two Tubes

To illustrate that the CPU-time spent on contact searching in relation to the total time spent on the solution of the structural equations does not increase with the increase of the number of nodes in the contact system, a problem with two tubes impacting each other is analyzed with two different mesh densities models. The diameter of each of the tubes is 0.20 m, thickness is 0.005 m, length is 0.46 m. The material of the two tubes is assumed to be elastoplastic with a linear isotropic hardening. One tube has an initial velocity of 35 m/s, another has the same initial velocity in the opposite direction. In the first model, each of the tubes is modelled by 300 4-node shell elements. The deformed configuration at time 6ms is shown in Fig.7(a). Total CPU time used is 5642s. The contact searching algorithm requires 649s, about 11.5% of the total CPU time. In the second model, each of the tubes is modelled by 912 shell elements. The deformed configuration at time 6 ms is shown in Fig.7(b). Total CPU time used is 16385 s. The contact searching algorithm takes 11.7% of the total CPU time, i.e. 1917s.

Vol.16, No.4

Wang Fujun et al.: A Contact Searching Algorithm

381

(a) for first model (b) for second model Fig.7 The deformed configurations at time 6 ms 4.3 I m p a c t o f V e h i c l e F r a m e A g a i n s t a R i g i d W a l l The third example concerns the impact between a vehicle frame and a rigid wall[13]. The finite element model used for the vehicle frame is shown in Fig.8. The vehicle frame has an initial velocity of 4 8 k m / h . A lumped mass of 550kg is distributed on the second transversal beam to simulate the engine loading. The cell structure for this problem is defined such that only the frontal part of the vehicle frame is included in the contact searching process. The deformed configuration of the frontal part of the vehicle frame at time 18.8ms is shown in Fig.9. Figure 10 gives the permanent deformation of the frame end point

Fig.8 The finite element model of a vehicle frame

Fig.9 The deformed configuration of frontal part

Fig.10 The displacement compared to experiment [~a] result

382

ACTA MECHANICA SINICA (English Series)

2000

A (indicated in Fig.8) in x direction as compmed to the experimental results [l:~,. T t > proportion of time spent on contact searching to that on the total solution is 7.3%, while it is 13.9% for the conventional method that has the order of O(Nlog2N ). 5 CONCLUSIONS The searching algorithm presented in this paper deals with the global contact searching in the dynamic finite element analysis, using tile cell structure and the node linked-list. \~,~ have used this algorithm in solving different contact-impact problems with very good results that demonstrate that the algorithm is excellent for solving 3-dimensional contact-impact problems. In summary, the new contact searching algorithm enjoys the following benefits. (1) The problem of sorting and searching in three dimensions is transformed to a retrieving process from two one-dimensional arrays. This is particularly suitable for parallel computing. (2) The algorithm does efficient sem'ching. Its cost of contact searching is of the order of O(N), as compared with OINlog2N ) for other old algorithms. (3) It can handle a contact system with an arbitrary mesh layout so that it allows treatment of contact and impact conditions anywhere in one model, including both twosurface contact and self-contact, without needing a surface definition by the user. (4) It is rather simple and its implementation in any finite element program is straight forward. In fact, this algorithm has been utilized with fully automatic treatment of contact in explicit finite element codes.

REFERENCES 10ldenburg M, Nilsson L. The position code algorithm for contact searching. Int J Numer Meth Engng, 1994, 37:359~386 2 Benson D J, Hallquist JO. A single surface contact algorithm for the post-buckling analysis of shell structures. Comput Methods Appl Mech Engrg, 1990, 78:141--q63 3 Wang SP, Nakamachi E. The inside-outside contact search algorithm for finite element analysis. Int J Numer Meth En.gng, 1997, 40(19): 3665~3685 4 Hallquist JO. LS-DYNA3D Theoretical Manual. Livermore: Livermore Software Technology Corporation, 1993. 135~180 5 Hallquist JO, Goudreau GL, Bension DJ. Sliding interfaces with contact-impact in large-scale Lagrangian computations. Comput Methods Appl Meeh Engrg, 1985, 51:107,-o137 6 Munjiza A, Andrews KRF. NBS contact detection algorithm for bodies of similar size. Int J Numer Meth Engng, 1998, 43:131~149 7 Yah WM. Data Structure. Beijing: Tsinghua University Press, 1992. 10~34 (in Chinese) 8 Bathe KJ. Finite Element Procedures. Englewood Cliffs: Prentice-Hall, 1996. 98~120 9 Zhong ZH. Finite Element Procedures for Contact-Impact Problems. Oxford: Oxford University Press, 1993. 287~313 10 Belytschko T, Neal MO. Contact-impact by the pinball algorithm with penalty and Langrangian methods. Int J Numer Meth Engng, 1991, 31:547~572 11 Belytschko T, Lin JI. A three-dimensional impact-penetration algorithm with erosion. Computers and Structures, 1987, 25(1): 95~104 12 Bathe K J, Walczak JW, Guillermin O, et al. Advances in crush analysis. Computers and Structures, 1999, 72(1~3): 31~47 13 \Vang CY. Study on the deformation and dynamic response of vehicle frame in crash. [dissertation]. Beijing: Tsinghua University, 1996. 55~60