- Email: [email protected]
A modified tree search algorithm for contact BEM-LCP approach
Pergamon
Computers & S!rucrures Vol. 59, No. 4, pp. 707-713, 1996 Copyright CI 1996 Elsevw Science Ltd Printed in Great Britain. All rights reserved 0045-7949196 $I 5.00 + 0.00
00457949(95)00286-3
A MODIFIED TREE SEARCH ALGORITHM CONTACT BEM-LCP APPROACH
FOR
B. Q. Xu and X. A. Kong Institute
of Applied
Mechanics,
Southwest Jiaotong University, People’s Republic of China (Received
30 September
Chengdu,
Sichuan,
610031,
1994)
Abstract-The present article is mainly concerned with the modification of the tree search algorithm for the solution of BEM-LCP (boundary elements method-linear complementary programming) equation emerged in general frictional contact problems. Numerical procedure, herein indicated, features a suggestion that the normal and tangential tractions are considered as a whole system, such that the
influence of tangential on the normal contact variables can be included. Modification and innovation on search path are made in the present work to alleviate the storage difficulties encountered in existing tree search algorithm. Numerical examinations further justifies the effectiveness of the present algorithm, in specific, in the solution of “non-quasi-identical” contact problems.
brief discussion is provided by Klarbring [8]. According to Cottle and Jahanshahlou, the modified simplex method is generally valid for positive definite matrices, principal pivoting method (p.p.m) plausible to P-matrix, modified p.p.m to positive semi-definite matrix, and, finally, the Lemke method is suitable to co-positive plus matrices, whereas the matrix related to LCP deduced in a BEM-LCP scheme is non-symmetric and positive indefinite [ 141. Its structure is not quite clear and is still to be investigated, according to published literature. Facts demonstrated that the widely used Lemke algorithm is not quite suitable to the BEM-PLCP matrix: either it does not provide access to the solution, or the iteration process is unstable. In order to overcome this numerical difficulty, Kong et al. [15] developed a two stage algorithm by separating the normal from the tangential contact parts and having them treated independently. However, this technique is basically effective for “quasi-identical” contact problems [16], as it is not practical to have it extended to situations where the influence of tangential contact tractions on the normal contact pressures is not negligible. The present work is aimed at establishing a solution algorithm suitable to a non-symmetric, positive indefinite matrix related to the BEM-LCP contact analysis which is basically not quasi-identical. For this purpose, the tree search algorithm proposed in Ref. [l l] is invoked. Attractive in generality, as indicated by its originators, in as much as being free from any assumptions on the characteristics of the LCP matrix and of providing every solution of the LCP, its adaptability to engineering applications needs further investigation, so far as basic research in mathematical programming is concerned. Even worse, the algorithm is fatally disadvantaged in the
1. INTRODUCTION
The strong non-linearity exhibited in contact problems stimulates the invocation of two kinds of solution procedures: the trial-and-error algorithms and algorithms related to mathematical programming (MP) formulations. As indicated by some authors, the trial-and-error algorithm cannot always ensure to capture the non-unique nature of frictional contact problems and may fail to converge [l]. On the other hand, frictional contact problems, either with constant or varying contact surfaces, can be formulated as a unified LCP scheme of mathematical programming, once the three-dimensional Coulomb’s friction conical surface is approximated by a polytope. According to recent publications, this numerical scheme is mainly a combination of linear complementarity formulation of elastic (or elasto-plastic) contact problems with the FE discretization [2-91; and the resulted FEM-LCPs are mostly solved by the principal pivoting [5-81, modified simplex [2, 31, and Lemke method [9]. Among them, the most popular methods for the solution of LCPs are summarized and discussed by Cottle [lo] and Jahanshahlou [ 111. Recently, a standard direct BEM formulation coupled with mathematical programming is employed by Kwak [ 121 for analysis of two-dimensional frictional contact problems where the modified simplex method is invoked for the resulted LCP. Extension of the BEM-LCP method to more general three-dimensional situations is achieved by Gakwaya [l] and Kwak [13]. Apparently the advantages of the BEM based contact analysis over a FEM counterparts are that the effectiveness of algorithms for the solution the LCP arisen from a BEM discretization is still dubious. On the solvability of FEM-type contact LCP formulations, a 707
708
B. Q. XII and X. A. Kong
enormous memory space required, which makes it nearly impossible to have the algorithm implemented. As a matter of fact, it is required in the original algorithm to have all simplex tableaux stored, a fact that unavoidably demands unpredictably enormous memory space. Numerical experiences indicate that the number of simplex tableaux generated in the search process increase with an identical function of the number of variables. It is evident that without drastic modifications, application of existing tree search algorithm to engineering problems, such as contact problem, is practically impossible. In view of this, the tree search algorithm is immensely and drastically modified in the present study. Section 2 presents a brief formulation for the general contact problem of two elastic bodies with the introduction of total strain method. The modified tree search process is presented in Section 3. Emphasis is placed on the improvement of search path, such that the amazingly enormous computer memory requirement of the original algorithm can be brought down. To justify and assess the effectiveness of the present algorithm, numerical examinations are also conducted. It turns out that the present results are in good agreement with those of existing counterparts. 2.
BEM-LCP
FORMULATIONS FOR ELASTIC PROBLEMS
CONTACT
In the subsequent context, a brief description of the BEM-PLCP approach to elastic contact problems is presented. For detailed formulations the readers are referred to Ref. [l]. The standard direct boundary element formulation involving contact conditions and global equilibrium conditions leads to the discretized form of linear BEM system for each body. The unknowns outside the potential contact zone can be eliminated by a flexibility process. Fundamental equations for boundary element contact analysis are obtained with the law of action and reaction. On the potential contact boundary S,, the classical unilateral contact condition and Coulomb’s friction law are assumed to hold. For the case where the rigid body displacements are specified, according to Ref. [l]. the LCP can be described as is)
= {G) + {R)(t)
+ {MH*)
@NJ > 0, {PN} 3 0, {DN}T(P,V} = 0 {A} > 0, { 4)
> 0, {n}T{@} = 0,
contact problem can be solved by total stain method which is employed in the present paper. 3. MODIFIED
TREE SEARCH
METHOD
FOR LCP
According to Ref. [15], the non-standard LCP deduced in the preceding section cannot be effectively handled by the widely accepted methods such as principal pivoting, the Lemke method and the sort. Jahanshahlou and Mitra proposed a tree search method to deal with general LCP problems. The method is famed for its attractive generality in as much as no assumption on the characteristics of the matrix is involved. Computer implementation of the tree search algorithm exposes practically the inapplicability of this general purpose method to contact problem, even if the dimension of system matrix is very limited. This is attributed to the huge amount of computer memory space required, as in the aforementioned method, all the simplex tableaux in the search procedure have to be kept, and the number of simplex tableaux is proportional to the factorial of the number of complementary variable pairs. It is therefore almost impossible to apply the tree search algorithm proposed by Jahanshahlou, without serious modification, to contact analysis and other engineering problems. The main modifications made to the tree search algorithm are listed below. 3.1. Data structure and search procedure from leaf to root In the original tree search algorithm, every simplex tableau must be memorized in the solution process. The large number of variables involved in BEM-LCP in contact problems renders enormous computer memory to be occupied and the required memory space unpredictable. In view of this, in the present work, only the index matrix (denoted by MBA in the subsequent context) of the simplex tableau and the current tableau in search process is recorded, the preceding tableau being recovered through an inverse transformation of pivot procedure. The integer
(1) (2) (3)
where {S} = {DN, -@}‘, {A} = {PN, A)? Explicit forms of matrix [G], [R] and [M] are available in Ref. [l], where they are constant matrices. In case external loading varies proportionally, {t} can be expressed as a constant vector multiplied by a varying scalar load factor, so the
Fig. 1. Search path in Jahanshahlou’s algorithm.
A modified tree search algorithm
709 Table 2 1
ErIEr p result
0.0 S
0.05 P
0.1 S
0.2 S
0.1
10
0.0 0.2 FFPF
0.0 0.2
F = failed; S = success; P = poor precision. computers can meet the memory demand and the memory space demand is predictable even. A com-
parison of the present modified memory space demand with that of the original tree search method, both applied to the example in part 4, is given in Table 1, where NC = number of potential contact node pairs. 3.2. Weight the variables
Fig. 2. Search path in present algorithm.
matrix MBA is defined as: MBA(i, 1) record the non-basic variable index if MBA(i, 1) > 0, the nonbasic variable is s; MBA(i, 1) < 0 implies that A is non-basic variable. MBA(i, 2) is the basic variable index similar to MBA(i, 1). MBA(i, 3) and MBA(i, 4) record if the ith line is starred and ith column is flagged, respectively. As such, this record scheme effectively brings down the memory requirement of every tableau from n2 real numbers and 4n integers to only 4n integers, where n is the number of complementary variable pairs. As long as the search path of Jahanshahlou is a to and fro path from root to leaf (Fig. 1), it is always necessary to have every tableau generated in the search process recorded. Since the number of tableaux is around Cr”(n!)2,t( being the proportion of the positive elements of matrix M, the demand in memory space is disastrous. A remedial measure to come out of this dilemma is the use of a from-leaf-toroot one-way search path, as shown in Fig. 2. Note that in Figs 1 and 2, arrows are provided to indicate the search path. In fact, in this proposed search procedure, only one branch has to be recorded while the length of branch is about n. As a sharp contrast to the original algorithm, the number of tableaux that have to be recorded in the present situation is reduced from c~“(n!)~ to n. As a matter of fact, personal Table 1 Memory space occupied
NC
P
3
0.0
5 7
0.1 0.2 0.0 0.0
Jahanshahlou 9900 4950 14355 2.65 x lo5 6.25 x 10sa
Present 162 162 162 450 882
“The problem cannot be solved by the original tree search method, the demanded memory space is estimated through the search path of our modified algorithm.
In eqn (I), variables S,, A, are of different orders magnitude. Generally speaking, P,,, N @, PN/E <
m,;=m;E gj =g:E,
(i,j=l,NNTC)
ri = {Cr,-t,}*E
rnb = m,/E
(i, j = NNTC
(4)
(i = 1, NNTC)
(5)
+ 1, ZDZM)
(6)
{S}’ = {G}’ + {R}’ + [M]‘{A}’
(7)
{S}’ 2 0, {A}’ > 0, {S}‘T{A}’ = 0
(8)
DN=D;/E,I
=n’lE,
(9)
where NNTC is the number of potential contact node pairs. ZDZM is the number of complementary variable pairs. E is Young’s modulus of one of the contact bodies. In the search procedure, elements of M are set to 0 if their absolute values are very small, say, less than lo-‘, to get rid of the accumulation of computational error. 3.3. Accelerate search procedure In the original tree search algorithm, the column selection is carried out either by eliminating the basic complementary variable pair, or picking up the column with the minimum column index number. In as much as all the columns are of equal opportunity in the selection, the iteration procedure in searching for final solutions is unavoidably elongated and thus time consuming.
710
B. Q. Xu and X. A. Kong
I Weight
variable
L=O, L 1=L+
START
Basic feasible
1 Tab-O
Select column
4
AUX
F
Y
L=L-1
Pivot S
I
Ll=L+I
Reserve
Tab-L
J
Fig. 3. Flow chart of Jahanshahlou’s S = success.
algorithm.
F = failure;
Fig.
4.
Flow
chart
of present algorithm. S = success.
F = failure;
A modified tree search algorithm
711
3 ty ,
1OOmm 1111 1111
q=l.2N/
7
mm2
I;;i ‘2 0.5
h
9 +
0.3
F 3 d z
0.0 0.0
B
A
a ..+c
0.3 Distance
1.0
0.8
0.5 from
centre
(x/50)
Fig. 5. Structure geometry and BE discretization.
Fig. 7. W, distribution.
Unlike the original method, the auxiliary function is assumed in the present work as
Step 0. Weight the variables, i.e. to transform eqn (1) into eqn (7). Step 1. Apply phase I of the simplex method to system (7) to find a basic feasible solution. If there is no basic feasible solution of eqn (7), then there is no solution for the fundamental problem, and go to step 4; otherwise number the tableau associated with the basic feasible solution as Tableau-O, set L = 0. Record basic variables and non-basic variables in MBA (0). Step 2. Pick Tableau-L and choose the column, as described in Section 3.3. If it fails, go to auxiliary sequence, otherwise go to pivot procedure. If a pivotal transformation is carried out, set L = L + 1. Record basic variables and non-basic variables in MBA (L + 1) and go to step 2. If the auxiliary sequence ends in the terminal step, go to step 3.
f=lSl.lAl-
1
si- 1
s, < 0
(10)
Ai.
Ai
It is obvious thatf > 0 and with the unique minimum 0. Since the solution off= 0 is equivalent to the solution of the LCP [eqns (7H9)], the search process can be accelerated by selecting the column which can decrease the value of f with as great a speed as possible. 3.4. The modified tree search algorithm The steps of the aforementioned stated as follows:
algorithm can be
2.0
p&q
-Takahash -Kong
0.0
0.3
Distance
0.5 ‘from
0.8 centre
Fig. 6. PN and P, distribution.
1 .b (x/50)
distance
from
centre
(x/50)
Fig. 8. P, distribution for different values of Ep/EF.
712
B. Q. Xu and X. A. Kong
Step 3. If L=O go to step 4. Set L=L-l. Reserve tableau L and go to step 2. Step 4. The tree search is completed. The auxiliary sequence can be found in Jahanshahlou algorithm [ 111. Flow charts of the original and present modified tree search algorithm are illustrated in Figs 3 and 4, respectively. 4. NUMERICAL
EXAMPLE
A Fortran code module is incorporated into a general purpose BEM program BEFC running on the Siemens-7570.e computer environment. A numerical example is provided to justify the efficiency of the proposed algorithm. The extent of memory space saved due to our modification is investigated. It turns out that the attainment of dramatic savings in the amount of memory space is achieved via the invocation of the present modified algorithm. The influence of material properties of contact bodies is also investigated. Example. Contact between flat surfaces with different material properties In this classical two-dimensional conforming contact problem, the memory space saving capability of our modified tree search algorithm, as compared with the original tree search method, is studied. The interest is also in investigating the ability of the algorithm to deal with non quasi-identical contact problem [16]. The indented punch and foundation are illustrated in Fig. 5. The material constants are: Poisson’s ratio v = 0.35 for the two bodies, and three combinations of EP (elastic modulus of punch) and E, (elastic modulus of foundation) are taken as follows:
EP = 4000 N mm-‘,
E, = 40,000 N mm-’
Table 1 compares the corresponding memory space occupied by the original algorithm against the present one. It can be seen that our modification decreases the demanded memory space from a”@ !)’ to 2n 2. This is much less than that demanded by the Jahanshahlou algorithm. Another feature is that the required memory volume can be predicted in our algorithm, while, it is not predictable in the original one. Table 2 reveals that the solution precision cannot always be guaranteed without the present modification. Figure 6 shows the contact pressure and tangential distribution for p = 0.2, EP/EF = 1. The results are compared with the flexibility approach solutions of Ref. [17] and two stage solutions [14]. The present solution is in good agreement with that of Kong. The
singularity at the contact edge, which appears in Takahash solutions, no longer exists. The relative tangential displacement distributions for various friction coefficients are shown in Fig. 7. The present solution varies more smoothly on the contact margin. Finally, Fig. 8 gives the distribution of P, for different material combinations (p = 0.2). It can be seen that high ratio of EP/EF makes larger sticking region, which is also observed independently by Takahash and by Kong. The negative tangential tractions are also computed for E,/E, = 10. The negative tangential displacements obtained by Takahash, disadvantaged in physical background, do not exist in our results. In the present example, all points in the contact zone are sticking as EP/EF = 10. 5. SUMMARY
The tree search algorithm is greatly improved for normal and tangential coupling contact BEM-LCP formulation in the present article. The emphasis is placed on overcoming the huge memory difficulty. The present modified tree search algorithm is valid for more general problems, such as frictional, nonquasi-identical contact problems, than the two stage algorithm used by Kong [ 151. The main features of the present algorithm include: (1) The tree search algorithm is modified and implemented into a BEM-LCP approach to general contact problem for the solution of the non-symmetric, positive indefinite LCP equation. The huge memory difficulty encountered in the original version of the tree search algorithm no longer exists in the present algorithm, and the solution precision is greatly improved. (2) The normal and tangential coupling contact problem can be solved simultaneously, which captures the real physical feature of contact problems. (3) By piecewise linear approximation of Coulomb’s friction law, the contact problem can be solved by the total strain method, if the external loads vary proportionally. In this way, the contact status iteration is no longer needed. (4) The so-called “non quasi-identical” contact problem can be handled in the same way as the “quasi-identical” problem. So, the influence of material property difference of contact bodies can be guaranteed. Acknowledgements-The authors would like to acknowledge the financial support from the Foundation of Education Commission of China. The authors also wish to express their gratitude to Professor Chen for his helpful comments. REFERENCES
1. A. Gakwaya, D. Lambert element and mathematical
and A. Cardou, A boundary programming approach for
A modified
2.
3.
4. 5.
6.
7.
8.
9.
tree search
frictional contact problems. Comput. Struct. 42, 341-353 (1992). R. Chand, E. J. Haug and K. Rim, Analysis of unbonded contact problems by means of quadratic programming. J. Optim. Ther. Appl. 20, 171-189 (1976). T. F. Conry and A. Seireg, A mathematical programming method for design of elastic bodies in contact. J. appl. Mech. 387-392 (1971). J. J. Kalker, Contact mechanical algorithms. Commun. appl. numer. Meth. 4, 25-32 (1984). G. Holmberg, A solution scheme for three-dimensional multi-body contact problem with friction using mathematical programming. Comput. Struct. 37, 503-514 (1988). A. Klarbring, General contact boundary conditions and the analysis of frictional systems. Int. J. Solidr Struct. 22, 1377-1398 (1986). A. Klarbring, A mathematical programming approach to three-dimensional contact problem with friction. Comput. Meth. appl. Mech. Engng 58, 175-200 (1986). A. Klarbring and G. Bjorkman, A mathematical programming approach to contact problem with friction and varying contact surface. Comput. Struct. 30, 1185-1198 (1988). S. M. Sun and W. X. Zhong, A finite element method for elasto-plastic structure and contact problem by parametric quadratic programming. Znt. J. numer. Meth. Engng 26, 2723-2738 (1989).
algorithm
713
pivot 10. R. W. Cottle and G. B. Dantzig, Complementary theory of mathematic programming. Linear Algebra A&. 1, 103-125 (1968). and G. Mitra, Linear complemen11. G. R. Jahanshahlou tary problem and a tree search algorithm for its solution. In: Proc. 9th Ini. Mathematical Programming Symp., Budapest (1976). problem 12. B. M. Kwak and S. S. Lee. A comulementaritv formulation for two-dimensional frictional -contact problems. Comput. Sfruct. u), 469-480 (1988). problem formulation 13. B. M. Kwak, Complementarity for three-dimensional frictional contact. J. appl. Mech. 58, 134140 (1991). supplement of FEM. 14. Z. H. Yao, BEM is an important Lecture on Int. Modern Solid Mechanics Problems Study Groun. Beiing (1992). A. Cardou and L. Cloutier, 15. X. A. Kong; A. Gakwaya, A numerical solution of general frictional contact problems by the direct boundary element and mathematical programming approach. Comput. Strucr. 45, 95-l 12 (1993). 16. J. J. Kalker, Three-dimensional Elastic Bodies in Rolling Contact. Klumer Academic (1990). and C. A. Brebbia, Analysis of contact 17. S. Takahash problems in elastic bodies using a B.E.M. flexibility approach. In: Boundary Element Method X (Edited by C. A. Brebbia), pp. 353-379. Springer, Berlin (1988).
Computers & S!rucrures Vol. 59, No. 4, pp. 707-713, 1996 Copyright CI 1996 Elsevw Science Ltd Printed in Great Britain. All rights reserved 0045-7949196 $I 5.00 + 0.00
00457949(95)00286-3
A MODIFIED TREE SEARCH ALGORITHM CONTACT BEM-LCP APPROACH
FOR
B. Q. Xu and X. A. Kong Institute
of Applied
Mechanics,
Southwest Jiaotong University, People’s Republic of China (Received
30 September
Chengdu,
Sichuan,
610031,
1994)
Abstract-The present article is mainly concerned with the modification of the tree search algorithm for the solution of BEM-LCP (boundary elements method-linear complementary programming) equation emerged in general frictional contact problems. Numerical procedure, herein indicated, features a suggestion that the normal and tangential tractions are considered as a whole system, such that the
influence of tangential on the normal contact variables can be included. Modification and innovation on search path are made in the present work to alleviate the storage difficulties encountered in existing tree search algorithm. Numerical examinations further justifies the effectiveness of the present algorithm, in specific, in the solution of “non-quasi-identical” contact problems.
brief discussion is provided by Klarbring [8]. According to Cottle and Jahanshahlou, the modified simplex method is generally valid for positive definite matrices, principal pivoting method (p.p.m) plausible to P-matrix, modified p.p.m to positive semi-definite matrix, and, finally, the Lemke method is suitable to co-positive plus matrices, whereas the matrix related to LCP deduced in a BEM-LCP scheme is non-symmetric and positive indefinite [ 141. Its structure is not quite clear and is still to be investigated, according to published literature. Facts demonstrated that the widely used Lemke algorithm is not quite suitable to the BEM-PLCP matrix: either it does not provide access to the solution, or the iteration process is unstable. In order to overcome this numerical difficulty, Kong et al. [15] developed a two stage algorithm by separating the normal from the tangential contact parts and having them treated independently. However, this technique is basically effective for “quasi-identical” contact problems [16], as it is not practical to have it extended to situations where the influence of tangential contact tractions on the normal contact pressures is not negligible. The present work is aimed at establishing a solution algorithm suitable to a non-symmetric, positive indefinite matrix related to the BEM-LCP contact analysis which is basically not quasi-identical. For this purpose, the tree search algorithm proposed in Ref. [l l] is invoked. Attractive in generality, as indicated by its originators, in as much as being free from any assumptions on the characteristics of the LCP matrix and of providing every solution of the LCP, its adaptability to engineering applications needs further investigation, so far as basic research in mathematical programming is concerned. Even worse, the algorithm is fatally disadvantaged in the
1. INTRODUCTION
The strong non-linearity exhibited in contact problems stimulates the invocation of two kinds of solution procedures: the trial-and-error algorithms and algorithms related to mathematical programming (MP) formulations. As indicated by some authors, the trial-and-error algorithm cannot always ensure to capture the non-unique nature of frictional contact problems and may fail to converge [l]. On the other hand, frictional contact problems, either with constant or varying contact surfaces, can be formulated as a unified LCP scheme of mathematical programming, once the three-dimensional Coulomb’s friction conical surface is approximated by a polytope. According to recent publications, this numerical scheme is mainly a combination of linear complementarity formulation of elastic (or elasto-plastic) contact problems with the FE discretization [2-91; and the resulted FEM-LCPs are mostly solved by the principal pivoting [5-81, modified simplex [2, 31, and Lemke method [9]. Among them, the most popular methods for the solution of LCPs are summarized and discussed by Cottle [lo] and Jahanshahlou [ 111. Recently, a standard direct BEM formulation coupled with mathematical programming is employed by Kwak [ 121 for analysis of two-dimensional frictional contact problems where the modified simplex method is invoked for the resulted LCP. Extension of the BEM-LCP method to more general three-dimensional situations is achieved by Gakwaya [l] and Kwak [13]. Apparently the advantages of the BEM based contact analysis over a FEM counterparts are that the effectiveness of algorithms for the solution the LCP arisen from a BEM discretization is still dubious. On the solvability of FEM-type contact LCP formulations, a 707
708
B. Q. XII and X. A. Kong
enormous memory space required, which makes it nearly impossible to have the algorithm implemented. As a matter of fact, it is required in the original algorithm to have all simplex tableaux stored, a fact that unavoidably demands unpredictably enormous memory space. Numerical experiences indicate that the number of simplex tableaux generated in the search process increase with an identical function of the number of variables. It is evident that without drastic modifications, application of existing tree search algorithm to engineering problems, such as contact problem, is practically impossible. In view of this, the tree search algorithm is immensely and drastically modified in the present study. Section 2 presents a brief formulation for the general contact problem of two elastic bodies with the introduction of total strain method. The modified tree search process is presented in Section 3. Emphasis is placed on the improvement of search path, such that the amazingly enormous computer memory requirement of the original algorithm can be brought down. To justify and assess the effectiveness of the present algorithm, numerical examinations are also conducted. It turns out that the present results are in good agreement with those of existing counterparts. 2.
BEM-LCP
FORMULATIONS FOR ELASTIC PROBLEMS
CONTACT
In the subsequent context, a brief description of the BEM-PLCP approach to elastic contact problems is presented. For detailed formulations the readers are referred to Ref. [l]. The standard direct boundary element formulation involving contact conditions and global equilibrium conditions leads to the discretized form of linear BEM system for each body. The unknowns outside the potential contact zone can be eliminated by a flexibility process. Fundamental equations for boundary element contact analysis are obtained with the law of action and reaction. On the potential contact boundary S,, the classical unilateral contact condition and Coulomb’s friction law are assumed to hold. For the case where the rigid body displacements are specified, according to Ref. [l]. the LCP can be described as is)
= {G) + {R)(t)
+ {MH*)
@NJ > 0, {PN} 3 0, {DN}T(P,V} = 0 {A} > 0, { 4)
> 0, {n}T{@} = 0,
contact problem can be solved by total stain method which is employed in the present paper. 3. MODIFIED
TREE SEARCH
METHOD
FOR LCP
According to Ref. [15], the non-standard LCP deduced in the preceding section cannot be effectively handled by the widely accepted methods such as principal pivoting, the Lemke method and the sort. Jahanshahlou and Mitra proposed a tree search method to deal with general LCP problems. The method is famed for its attractive generality in as much as no assumption on the characteristics of the matrix is involved. Computer implementation of the tree search algorithm exposes practically the inapplicability of this general purpose method to contact problem, even if the dimension of system matrix is very limited. This is attributed to the huge amount of computer memory space required, as in the aforementioned method, all the simplex tableaux in the search procedure have to be kept, and the number of simplex tableaux is proportional to the factorial of the number of complementary variable pairs. It is therefore almost impossible to apply the tree search algorithm proposed by Jahanshahlou, without serious modification, to contact analysis and other engineering problems. The main modifications made to the tree search algorithm are listed below. 3.1. Data structure and search procedure from leaf to root In the original tree search algorithm, every simplex tableau must be memorized in the solution process. The large number of variables involved in BEM-LCP in contact problems renders enormous computer memory to be occupied and the required memory space unpredictable. In view of this, in the present work, only the index matrix (denoted by MBA in the subsequent context) of the simplex tableau and the current tableau in search process is recorded, the preceding tableau being recovered through an inverse transformation of pivot procedure. The integer
(1) (2) (3)
where {S} = {DN, -@}‘, {A} = {PN, A)? Explicit forms of matrix [G], [R] and [M] are available in Ref. [l], where they are constant matrices. In case external loading varies proportionally, {t} can be expressed as a constant vector multiplied by a varying scalar load factor, so the
Fig. 1. Search path in Jahanshahlou’s algorithm.
A modified tree search algorithm
709 Table 2 1
ErIEr p result
0.0 S
0.05 P
0.1 S
0.2 S
0.1
10
0.0 0.2 FFPF
0.0 0.2
F = failed; S = success; P = poor precision. computers can meet the memory demand and the memory space demand is predictable even. A com-
parison of the present modified memory space demand with that of the original tree search method, both applied to the example in part 4, is given in Table 1, where NC = number of potential contact node pairs. 3.2. Weight the variables
Fig. 2. Search path in present algorithm.
matrix MBA is defined as: MBA(i, 1) record the non-basic variable index if MBA(i, 1) > 0, the nonbasic variable is s; MBA(i, 1) < 0 implies that A is non-basic variable. MBA(i, 2) is the basic variable index similar to MBA(i, 1). MBA(i, 3) and MBA(i, 4) record if the ith line is starred and ith column is flagged, respectively. As such, this record scheme effectively brings down the memory requirement of every tableau from n2 real numbers and 4n integers to only 4n integers, where n is the number of complementary variable pairs. As long as the search path of Jahanshahlou is a to and fro path from root to leaf (Fig. 1), it is always necessary to have every tableau generated in the search process recorded. Since the number of tableaux is around Cr”(n!)2,t( being the proportion of the positive elements of matrix M, the demand in memory space is disastrous. A remedial measure to come out of this dilemma is the use of a from-leaf-toroot one-way search path, as shown in Fig. 2. Note that in Figs 1 and 2, arrows are provided to indicate the search path. In fact, in this proposed search procedure, only one branch has to be recorded while the length of branch is about n. As a sharp contrast to the original algorithm, the number of tableaux that have to be recorded in the present situation is reduced from c~“(n!)~ to n. As a matter of fact, personal Table 1 Memory space occupied
NC
P
3
0.0
5 7
0.1 0.2 0.0 0.0
Jahanshahlou 9900 4950 14355 2.65 x lo5 6.25 x 10sa
Present 162 162 162 450 882
“The problem cannot be solved by the original tree search method, the demanded memory space is estimated through the search path of our modified algorithm.
In eqn (I), variables S,, A, are of different orders magnitude. Generally speaking, P,,, N @, PN/E <
m,;=m;E gj =g:E,
(i,j=l,NNTC)
ri = {Cr,-t,}*E
rnb = m,/E
(i, j = NNTC
(4)
(i = 1, NNTC)
(5)
+ 1, ZDZM)
(6)
{S}’ = {G}’ + {R}’ + [M]‘{A}’
(7)
{S}’ 2 0, {A}’ > 0, {S}‘T{A}’ = 0
(8)
DN=D;/E,I
=n’lE,
(9)
where NNTC is the number of potential contact node pairs. ZDZM is the number of complementary variable pairs. E is Young’s modulus of one of the contact bodies. In the search procedure, elements of M are set to 0 if their absolute values are very small, say, less than lo-‘, to get rid of the accumulation of computational error. 3.3. Accelerate search procedure In the original tree search algorithm, the column selection is carried out either by eliminating the basic complementary variable pair, or picking up the column with the minimum column index number. In as much as all the columns are of equal opportunity in the selection, the iteration procedure in searching for final solutions is unavoidably elongated and thus time consuming.
710
B. Q. Xu and X. A. Kong
I Weight
variable
L=O, L 1=L+
START
Basic feasible
1 Tab-O
Select column
4
AUX
F
Y
L=L-1
Pivot S
I
Ll=L+I
Reserve
Tab-L
J
Fig. 3. Flow chart of Jahanshahlou’s S = success.
algorithm.
F = failure;
Fig.
4.
Flow
chart
of present algorithm. S = success.
F = failure;
A modified tree search algorithm
711
3 ty ,
1OOmm 1111 1111
q=l.2N/
7
mm2
I;;i ‘2 0.5
h
9 +
0.3
F 3 d z
0.0 0.0
B
A
a ..+c
0.3 Distance
1.0
0.8
0.5 from
centre
(x/50)
Fig. 5. Structure geometry and BE discretization.
Fig. 7. W, distribution.
Unlike the original method, the auxiliary function is assumed in the present work as
Step 0. Weight the variables, i.e. to transform eqn (1) into eqn (7). Step 1. Apply phase I of the simplex method to system (7) to find a basic feasible solution. If there is no basic feasible solution of eqn (7), then there is no solution for the fundamental problem, and go to step 4; otherwise number the tableau associated with the basic feasible solution as Tableau-O, set L = 0. Record basic variables and non-basic variables in MBA (0). Step 2. Pick Tableau-L and choose the column, as described in Section 3.3. If it fails, go to auxiliary sequence, otherwise go to pivot procedure. If a pivotal transformation is carried out, set L = L + 1. Record basic variables and non-basic variables in MBA (L + 1) and go to step 2. If the auxiliary sequence ends in the terminal step, go to step 3.
f=lSl.lAl-
1
si- 1
s, < 0
(10)
Ai.
Ai
It is obvious thatf > 0 and with the unique minimum 0. Since the solution off= 0 is equivalent to the solution of the LCP [eqns (7H9)], the search process can be accelerated by selecting the column which can decrease the value of f with as great a speed as possible. 3.4. The modified tree search algorithm The steps of the aforementioned stated as follows:
algorithm can be
2.0
p&q
-Takahash -Kong
0.0
0.3
Distance
0.5 ‘from
0.8 centre
Fig. 6. PN and P, distribution.
1 .b (x/50)
distance
from
centre
(x/50)
Fig. 8. P, distribution for different values of Ep/EF.
712
B. Q. Xu and X. A. Kong
Step 3. If L=O go to step 4. Set L=L-l. Reserve tableau L and go to step 2. Step 4. The tree search is completed. The auxiliary sequence can be found in Jahanshahlou algorithm [ 111. Flow charts of the original and present modified tree search algorithm are illustrated in Figs 3 and 4, respectively. 4. NUMERICAL
EXAMPLE
A Fortran code module is incorporated into a general purpose BEM program BEFC running on the Siemens-7570.e computer environment. A numerical example is provided to justify the efficiency of the proposed algorithm. The extent of memory space saved due to our modification is investigated. It turns out that the attainment of dramatic savings in the amount of memory space is achieved via the invocation of the present modified algorithm. The influence of material properties of contact bodies is also investigated. Example. Contact between flat surfaces with different material properties In this classical two-dimensional conforming contact problem, the memory space saving capability of our modified tree search algorithm, as compared with the original tree search method, is studied. The interest is also in investigating the ability of the algorithm to deal with non quasi-identical contact problem [16]. The indented punch and foundation are illustrated in Fig. 5. The material constants are: Poisson’s ratio v = 0.35 for the two bodies, and three combinations of EP (elastic modulus of punch) and E, (elastic modulus of foundation) are taken as follows:
EP = 4000 N mm-‘,
E, = 40,000 N mm-’
Table 1 compares the corresponding memory space occupied by the original algorithm against the present one. It can be seen that our modification decreases the demanded memory space from a”@ !)’ to 2n 2. This is much less than that demanded by the Jahanshahlou algorithm. Another feature is that the required memory volume can be predicted in our algorithm, while, it is not predictable in the original one. Table 2 reveals that the solution precision cannot always be guaranteed without the present modification. Figure 6 shows the contact pressure and tangential distribution for p = 0.2, EP/EF = 1. The results are compared with the flexibility approach solutions of Ref. [17] and two stage solutions [14]. The present solution is in good agreement with that of Kong. The
singularity at the contact edge, which appears in Takahash solutions, no longer exists. The relative tangential displacement distributions for various friction coefficients are shown in Fig. 7. The present solution varies more smoothly on the contact margin. Finally, Fig. 8 gives the distribution of P, for different material combinations (p = 0.2). It can be seen that high ratio of EP/EF makes larger sticking region, which is also observed independently by Takahash and by Kong. The negative tangential tractions are also computed for E,/E, = 10. The negative tangential displacements obtained by Takahash, disadvantaged in physical background, do not exist in our results. In the present example, all points in the contact zone are sticking as EP/EF = 10. 5. SUMMARY
The tree search algorithm is greatly improved for normal and tangential coupling contact BEM-LCP formulation in the present article. The emphasis is placed on overcoming the huge memory difficulty. The present modified tree search algorithm is valid for more general problems, such as frictional, nonquasi-identical contact problems, than the two stage algorithm used by Kong [ 151. The main features of the present algorithm include: (1) The tree search algorithm is modified and implemented into a BEM-LCP approach to general contact problem for the solution of the non-symmetric, positive indefinite LCP equation. The huge memory difficulty encountered in the original version of the tree search algorithm no longer exists in the present algorithm, and the solution precision is greatly improved. (2) The normal and tangential coupling contact problem can be solved simultaneously, which captures the real physical feature of contact problems. (3) By piecewise linear approximation of Coulomb’s friction law, the contact problem can be solved by the total strain method, if the external loads vary proportionally. In this way, the contact status iteration is no longer needed. (4) The so-called “non quasi-identical” contact problem can be handled in the same way as the “quasi-identical” problem. So, the influence of material property difference of contact bodies can be guaranteed. Acknowledgements-The authors would like to acknowledge the financial support from the Foundation of Education Commission of China. The authors also wish to express their gratitude to Professor Chen for his helpful comments. REFERENCES
1. A. Gakwaya, D. Lambert element and mathematical
and A. Cardou, A boundary programming approach for
A modified
2.
3.
4. 5.
6.
7.
8.
9.
tree search
frictional contact problems. Comput. Struct. 42, 341-353 (1992). R. Chand, E. J. Haug and K. Rim, Analysis of unbonded contact problems by means of quadratic programming. J. Optim. Ther. Appl. 20, 171-189 (1976). T. F. Conry and A. Seireg, A mathematical programming method for design of elastic bodies in contact. J. appl. Mech. 387-392 (1971). J. J. Kalker, Contact mechanical algorithms. Commun. appl. numer. Meth. 4, 25-32 (1984). G. Holmberg, A solution scheme for three-dimensional multi-body contact problem with friction using mathematical programming. Comput. Struct. 37, 503-514 (1988). A. Klarbring, General contact boundary conditions and the analysis of frictional systems. Int. J. Solidr Struct. 22, 1377-1398 (1986). A. Klarbring, A mathematical programming approach to three-dimensional contact problem with friction. Comput. Meth. appl. Mech. Engng 58, 175-200 (1986). A. Klarbring and G. Bjorkman, A mathematical programming approach to contact problem with friction and varying contact surface. Comput. Struct. 30, 1185-1198 (1988). S. M. Sun and W. X. Zhong, A finite element method for elasto-plastic structure and contact problem by parametric quadratic programming. Znt. J. numer. Meth. Engng 26, 2723-2738 (1989).
algorithm
713
pivot 10. R. W. Cottle and G. B. Dantzig, Complementary theory of mathematic programming. Linear Algebra A&. 1, 103-125 (1968). and G. Mitra, Linear complemen11. G. R. Jahanshahlou tary problem and a tree search algorithm for its solution. In: Proc. 9th Ini. Mathematical Programming Symp., Budapest (1976). problem 12. B. M. Kwak and S. S. Lee. A comulementaritv formulation for two-dimensional frictional -contact problems. Comput. Sfruct. u), 469-480 (1988). problem formulation 13. B. M. Kwak, Complementarity for three-dimensional frictional contact. J. appl. Mech. 58, 134140 (1991). supplement of FEM. 14. Z. H. Yao, BEM is an important Lecture on Int. Modern Solid Mechanics Problems Study Groun. Beiing (1992). A. Cardou and L. Cloutier, 15. X. A. Kong; A. Gakwaya, A numerical solution of general frictional contact problems by the direct boundary element and mathematical programming approach. Comput. Strucr. 45, 95-l 12 (1993). 16. J. J. Kalker, Three-dimensional Elastic Bodies in Rolling Contact. Klumer Academic (1990). and C. A. Brebbia, Analysis of contact 17. S. Takahash problems in elastic bodies using a B.E.M. flexibility approach. In: Boundary Element Method X (Edited by C. A. Brebbia), pp. 353-379. Springer, Berlin (1988).