On the rigid body displacements and rotations in unilateral contact problems and applications

0045-7949191 $3.00 +

Computers & Strucrures Vol. 40, No. 3, pp. 599-514. 1991

Printedin GreatBritain.

0

0.00 1991 Pngamon Press plc

ON THE RIGID BODY DISPLACEMENTS AND ROTATIONS IN UNILATERAL CONTACT PROBLEMS AND APPLICATIONS G. E. STAVROULAKIS,~ P. D. PANAGIoToPouLost$ and A. M. AL-FAHED~ tSchoo1 of Technology, Aristotle University, GR-540 06 Thessaloniki, Greece $Institute of Technical Mechanics, R.W.T.H., D-5100 Aachen, F.R. Germany (Received 11 April 1990) Abstract-This paper deals with the treatment of rigid body displacements and rotations in unilateral contact problems. In the presence of rigid body modes the equivalent formulations of the problem, i.e. the variational inequality, the quadratic programming and the linear complementarity formulation involve positive semidefinite matrices. The basic differences from the classical bilateral problems is that in unilateral problems the rigid body modes must be compatible with the inequality constraints. This is a fairly difficult problem to solve and in its generality was open until now. In this paper a systematic analysis including convergence results is presented with respect to some methods already in use, namely the method of coupling with additional elastic bodies or ‘soft’ springs and the method of consideration of certain nodes as fixed. Moreover a new linear complementarity formulation of the problem which explicitly includes the rigid body ‘displacements’ is proposed and is numerically treated by a complementary pivoting technique. Necessary and sufficient conditions for the solution of the problem are derived and the theory is illustrated by examples from structural analysis and from robotics.

1. INTRODUCTION Problems involving inequalities in their constitutive relations or in their boundary conditions are frequently encountered in several engineering situations. This is the case of structures containing members capable of transmitting only certain type of stresses, e.g. tensionless or compressionless members, or subjected to noninterpenetration conditions along their boundaries or along prescribed interfaces as happens in tensionless supports in structures and crack-interfaces in rock masses. As well as the above problems yield criteria in the theories of plasticity and damage mechanics introduce inequalities in the complete formulation of the problem. All types of problems for which the principle of virtual work holds in inequality form, which is a variational inequality [l&4] or a hemivariational inequality [5-91 are called unilateral problems. Even in the range of small deformations and displacements and under linear elastic lxhaviour of the structure’s material, unilateral phenomena introduce high nonlinearities which cannot be dealt with satisfactorily by nonlinear structural analysis methods. To mention a few difficulties we recall that contact and noncontact regions in unilateral contact problems, plastic and elastic regions in plasticity, etc. are not known a priori and must be obtained from the numerical algorithm. Numerical and mathematical studies have shown that trial and error iterative techniques are time consuming and do not generally guarantee convergence to the correct

solution [3, 10-131. For this reason mathematical programming (MP) techniques have been applied in structural mechanics for the analysis of problems with inequalities [3, 14-161; they are accurate and quick. To give a measure of computational economy for 1000 contact unknowns the computational time of MP versus trial and error methods is l/20-1/25 [ 171. They are based on inequality constrained minimum propositions for the potential and/or the complementary energy which results from the variational formulation of the problem. In particular unilateral contact and/or friction boundary or interface relations for linearly elastic structures and in the framework of small deformations and displacements theory lead to quadratic programming problems (QPP) or equivalently to linear complementarity problems (LCP) [3, 12, 13,17-191. The same type of problems are obtained in discrete elastoplasticity according to Maier (see e.g. [20] for a complete reference list) and in the theory of cable structures [21] or of tensionless materials (cf. e.g. [3,22]). The subject of this paper is the derivation of several new results to complete the study of the rigid body displacements and rotations in unilateral contact problems and generally in all types of inequality problems which when discretized lead to a QPP or equivalently to an LCP. This type of inequality problems is most frequently encountered in several engineering applications and gives rise to the following QP problems [3, 18,23,24]

599

n(u) = min{x(v)/Av

- b < 0},

(1.1)

600

G.

E.

STAVROIJLAKIS et al.

where

2. METHODS FOR THE TREATMENT OF RIGID BODY DISPLACEMENTS AND ROTATIONS K

(II) = $uTKu - pTu

(1.2)

and n’(s) = min{n’(t)/Bt -c f 0},

(1.3)

n’(s) = fsTF,,s - s*e,.

(1.4)

where

Here a (xc) denotes the potential (complementary) energy, u, v (s, t) are the displacement (the stress) vectors, p is the loading-initial strain vector, K is the stiffness matrix, which is semidefinite if rigid body modes are not excluded, Au - b < 0 (Bt - c < 0) denotes the kinematically (statically) admissible set U,, (C,) for the displacements (the stresses), F, is the flexibility matrix and e, is the initial strain vector. Obviously (1.1) is equivalent (e.g. [25]) to the variational inequality uTK(v - u) - pT(v - u) 2 0,

Vv E U,,d.

(1.5)

Analogous equivalent variational inequality formulation exists for the problem (1.3). The treatment of rigid body displacements and rotations presents several difficulties even in classical bilateral linear problems because they give rise to positive semidefinite stiffness matrices. There are several approaches to the analysis of such problems under the presence of classical boundary conditions. However, in the presence of inequality subsidiary conditions, as happens, for example, in the unilateral contact problem with a rigid support, the problem of the study of such problems presents additional difficulties because the rigid body displacements or rotations must be compatible with the inequality constraints. The situation becomes more delicate in general unilateral problems expressed in terms of subdifferential boundary conditions of monotone or of nonmonotone type as happens, for example, in the case of ideally plastic hinges or in the case of plastic softening laws in frames and plates. All the aforementioned problems have been up to now dealt with by the following methods: (i) for a given loading we can add some fictitious springs of very low stiffness (‘soft springs’) and (ii) we can fix some nodes which seem to remain fixed during the loading process. All these approaches are systematized in the present paper, their validity is proved and they are compared with theoretical results obtained in the mathematical theory of inequality problems [2,3]. Moreover a new method is presented which permits the direct evaluation of the possible rigid body modes. This method is based on the complementary pivoting algorithm of Lemke. Finally the theory is illustrated by numerical results from structural analysis and from robotics.

2.1. Coupling with additiokal fictitious elastic bodies Let u E R* and let X be the set of solutions of the variational inequality problem (1.5). Since K is positive semidefinite the set X, assumed to be nonempty, is convex [26, p. 391 and closed for the Euclidean norm of the space B”. We denote each particular solution of (1.5) by u E X. Further a fictitious linear elastic structure is considered which contains the same set of discretization nodes and which has as a kinematically admissible set X. Moreover the structure has a positive definite stiffness matrix K, and is subjected to the loading ps. The equilibrium of this new structure is described by the variational inequality: Find II,,E X such as to satisfy the following variational inequality $K,(u*

- u,,) > p;(u* - II,),

Vu* E X.

(2.1)

Problem (2.1) has a unique solution II,,E X c V,,. The essential idea of the method proposed in this paragraph consists in the coupling of the initial structure with stiffness K with the additional fictitious elastic structure with stiffness KS by means of a real regularization parameter c > 0. Thus we are led to the study of a structure having a positive definite matrix K, = K + cK, which is subjected to the loading p6= p+ tps. The problem is to find II, E U,, such that ufK,(u* - u,) > p:(u* -II,),

Vu* E U,,.

(2.2)

Problem (2.2) has a unique solution due to the fact that Ks is positive definite. We shall show that the solution of problem (1.1) is the limit of the solution u, of (2.2) as L-+O+. We mention here that the idea of coupling an additional elastic structure has been extensively used, albeit in an intuitive way, in the calculation of highly nonlinear structural analysis phenomena, as, for example, in plastic analysis which takes into account softening and fracturing phenomena (compare in this respect the imbricate continuum model used by Baiant [27]). Concerning the behaviour of the model proposed the following proposition can be proved [26]. Proposition 2.1, Let X be nonempty, K be positive semidefinite and K, positive definite. We denote by II, the solution of (2.2) for each c > 0. Then as c-to+, u,+w, where w is a solution of the initial problem. Proof. Setting u* = II,,in (2.2) and u* = u,, u = u, in (1.5) and adding, we obtain

2 CP;(U~- n,),

(2.3)

Rigid body displacements in unilateral problems which, due to the fact that K is positive semidefinite and L is non-negative, gives ufK&

- u,) 2 P%, -uJ

uTK,u~%[[ull*,

V>O,

VuelR”,

where /.(j denotes the R” norm (Euclidean). (2.4) and (2.5) we obtain the relation p~u,+uTK,uo-p,T~,ZufK~u,~hPIIu,ll~

small positive number 6. Thus a positive definite stiffness matrix K, results K,=K+cI,

(2.4)

Since K, is positive definite we can determine number V > 0 such that

a

(2.5)

601

L >0

(2.13)

where I is the corresponding unit matrix. Let u be the solution of the unilateral problem (1.1). It is well known that the problem (1.1) is equivalent to the Karush-Kuhn-Tucker conditions [28] which give rise to the following nonstandard LCP

From

El-[-“*y][:]=[;‘]

(2.14)

(2.6) with

and thus from the Schwarz inequality we get Q II4 IVG Cl+

Y C2llU6

IL

where c, = -pib, c2 = IIK,u, + ps II are constants. Thus we have that

IIKII
20, r 2

0, y’r= 0.

(2.8)

where cj is a constant independent of c. Applying the Bolzano-Weirstrass theorem to (2.8) we may extract a subsequence, again denoted by {us} such that

Here y denotes the ‘gap’ of the unilateral constraint Av - b < 0 and r are the reactions generated by the inequality constraints (e.g. in a unilateral contact problem r are the reactions of the unilateral supports). The positive semidefinite matrix of the LCP (2.14) and (2.15) denoted by M will now be transformed into the positive definite matrix M,, M.=M+cI=[K_+;’

u,+w

in

Uad.

i.e. w E X.

Vu* E U,, ,

(2.10)

Q.E.D.

Remark. We can show further that w = u,, . Indeed

from (2.4) we get in the limit c +O+ wTKs(u, - w) 3 p;(u, -w>.

(2.11)

Since w o X setting u* = w in (2.1) and adding with (2.11) implies -(w-@‘K,(w-u,)>O

(2.12)

which implies that w = u,,. Proposition 2.1 shows that when the rigid body modes for a unilateral structure are suppressed by fictitious springs of very small stiffnesses, the obtained solution II, tends, as the stiffness of the additional springs tends to zero, to a solution of the initial structure. However, when additional springs are used, certain other approximation results can be derived. 2.2 Coupling with additional springs Let the stiffness of the fictitious springs be small (‘soft springs’); suppose that it is represented by a CAS 40,3--F

21.

(2.16)

(2.9)

Since U,, is closed, w E U,, . To complete the proof we must show that w is a solution of the initial problem, i.e. that w E X. Indeed taking the limit L + O+ in (2.2) we obtain that w’K(u* - w) > pT(u* - w),

(2.15)

(2.7)

Matrix M, corresponds to a structure having additional external springs to suppress the rigid body modes of the initial matrix K and to regularize the inequality constraints. For instance, in the unilateral contact problem with a rigid support, the additional springs give a Winkler-type compliance E to the rigid support, whereas in an ideal plasticity problem a small workhardening 6. The following propositions can be proved. Proposition 2.2. Suppose that the QPP (1.1) has a finite optimal solution. Then the LCP (2.14) and (2.15) corresponding to (1.1) admits a solution which is the limit, as c +O+, of the solution {II,, r,} of the 6 -perturbed LCP

[;]-$j=[ip]~ 3:

20,

rr 2 0,

yfr,=0.

(2.17) (2.18)

Proposition 2.3. Suppose that the QPP (1.1) has no solution, i.e. that either the inequality restrictions are infeasible or the potential energy of the structure is unbounded. Then as L +O+

These propositions are proved in the framework of quadratic programming theory in Ref. [29]. The reader is referred to this reference for the complete proof.

G. E.

602

STAVROULAKIS et

al.

2.3. Considering certain nodes as Jixed One commonly used technique for supporting parts of or whole structures which can exhibit rigid body displacements and rotations is to fix a sufficient number of nodes which we assume will fulfil the inequality constraints in equality form. For instance, in unilateral contact problems between deformable bodies one has to assume zero relative displacements of mutual pairs of interface nodes to prevent rigid body displacements or rotations of the whole ‘structure complex’ under consideration. In Fig. 1 one may specify that nodes K - K’ and A - A’ are tied, i.e. they do not debond. This is sufficient in order to fix the upper part of the whole structure and can be performed by using the well-known master-and-slave node technique. Then every available QPP program may be used in order to solve the unilateral contact problem. Additionally, one has to check after the solution if these nodes have been considered correctly. The simplest technique to do this is based on an extension for positive semidefinite matrices proved in Ref. [lo], of the third rule of the algorithm proposed by Theil and van de Panne for the solution of QP problems. This rule has been applied to the analysis of unilateral contact problems in Refs. [lo, 111. We denote by S the set of the inequality constraints which are fullfiled in equality form and let us be the solution of the classical bilateral problem n(uS) = min{x(v)/(Au -b), = 0,j E S}. The following proposition holds. Proposition 2.4. A solution us which satisfies all inequality restrictions of (1.1) coincides with the actual solution u, if and only if the solution ~‘-1~) violates the h-restriction, i.e. (Av - b), > 0 for every h ES. This proposition can be applied to check in our problem if the assumption that, e.g. K and K’ (Fig. 1) remain bonded is correct, after the calculation of the solution of the problem (see also [lo, 1I]). 3. ADDITIONAL SLACK VARIABLES AND LEMKE’S ALGORITHM FOR THE SOLUTION OF THE ARISING L.C.P. FORMULATION OF THE L.C.P.

In this section we shall derive a force-type method for the solution of the problem in which rigid body displacements are obtained as a byproduct of the calculation. We refer to the unilateral contact problem of Fig. 1 as a pilot problem where we assume, for the sake of simplicity, that the body B is subjected to classical (equality) boundary conditions. On B ‘rests’ another body A in which a node 0, hereafter called the pole of body A, is defined. The displacement vector of 0 completely describes the rigid body displacements and rotations of A. Note that any point of A can serve as the node 0. If the degrees of freedom used in a problem do not cover all rigid body

P..

Fig. 1. Schematic discretized problem, notation.

motions, additional degrees of freedom of other modes should be considered: for instance in plane elasticity problems having as degrees of freedom node displacement we must consider the two displacement degrees of freedom u, , u2of node 0 and the u2 degree of freedom of node 0, (cf. Fig. 1) in order to take into account the rigid body rotation of the body A. We introduce fictitious supports of A at the pole 0 (respectively at the poles 0, 0,) such as to give a statically determined substructure. Rigid body displacements and rotations are summed up to form a rigid body ‘displacement’ vector a, which for the three-dimensional problems reads

and analogously for other problems. The equilibrium equations for substructure A can now be written in the form Gr=Bp,.

(3.2)

Here G denotes the equilibrium matrix, which sums up the contributions of the unilateral nodal contact forces r, and B is a transformation matrix transforming the applied forces pA acting on A to the reference pole 0. Let us now consider the unilateral contact mechanism between parts A and B. Kinematic compatibility on the direction normal to the unilateral contact surface is written in the following form (cf. [30, 3 11) uk +

u; - 8; + b -

y = 0.

(3.3)

Here I$, uk denote the vectors of the contact displacements normal to the interface of the bodies A and 8, respectively and 8: denotes the normal displacements to the interface of the part A due to the rigid body ‘displacement’ d of 0. Vectors b and y are defined in the previous sections. After condensation we may write the elasticity laws uk=F Ar 7

I# = F,r,

(3.4)

for the parts A and B, respectively. Here F,, F, are the corresponding influence matrices, where F, is

603

Rigid body displacements in unilateral problems calculated for the isostatically supported substructure A at 0. At the interface we also assume that the tangential tractions have a prescribed value which is included in p*. From the kinematics of the part A which undergoes the rigid body displacement d at 0 we get the following expression -*=crft,

equal to zero in the final solution of the problem. Thus the following complementarity relations hold (v+)rI+ = 0,

(v-)%-

= 0.

(3.14)

Relations (3.7), (3.10a, b), (3.8), (3.11) and (3.14) may be written as the following LCP

(3.5) w=Mz+6,

Indeed the principle of complementary virtual work for part A gives for the real displacements f~, Ir, and for every statically admissible pair Bp,, r the following (BP,)~II + rTIi = 0.

w 20,

(3.15a)

z 2 0,

wrz = 0,

(3.15b)

where

(3.6)

Due to (3.2), we get from (3.6) rTGT&+ rTtIi = 0

for each r and thus (3.5) holds.

Equation (3.3) is written, by means of (3.4) and (3.5) in the form (F,+F,)r+GTI+b=y

(3.7)

and the unilateral contact law is written as in Section 2, i.e. r 2 0,

y 2 0,

rTy = 0.

(3.8)

Equations (3.2) and (3.7) with relations (3.8) form a nonstandard LCP which will be transformed below into a standard one by adding slack variables. To this end we first write equilibrium equations (3.2) as two inequalities GraBp,, -Gr

2 -BP,.

(3.9a) (3.9b)

By introducting additional nonnegative slack variables v+, v-, (3.9a, b) are written as equalities, i.e. Gr-v+ -Cr

=Bp,,

- v- = -BP,,

v+>o,

v-20.

(3.10a) (3.10b)

z;[:+.

b=[

:;I.

(3.1%)

To the standard LCP (3.15a-c) Lemke’s complementary pivoting, infeasibility reducing algorithm will be applied. Note at this point, that every unilateral problem of the type (1.1) or equivalently of the type (2.14) and (2. I5), may be written as a standard LCP like the one given by (3.15a-c) by splitting u into u+ and u-. However, this increases the number of unknowns in a nonacceptable way. Therefore the method of condensation as applied to the present pilot problem may be used, thus giving a standard LCP with a smaller number of unknowns. Actually we may always get an LCP with respect to a vector z including 8+, 8- and a number of unknowns r equal to the number of inequality restrictions. For instance, in the theory of piecewise linear yield surfaces in plasticity the dimension of r is equal to the number of finite elements (generally the number of Gauss integration points if finite elements of greater degree are used) multiplied by the number of hyperplanes linearizing the yield surface. Note that (3.15) can be obtained from the displacement method in a formalistic way: Kuhn-Tucker optimality conditions for QPP (I. 1) are written in the form

(3.11)

Ku+ATr=p,

(3.16a)

The rigid body ‘displacements’ is also written as a difference of two nonnegative vectors

Au+y=b,

(3.16b)

y>O, II=&‘-d-, fi+>o,

6-20.

r20,

yTr =O.

(3.17)

(3.12) (3.13)

Here the positive and the negative part of Q (i.e. 1;: = (6, + l&1)/2, J; = (-I& + [r&1)/2)is introduced. Slack variables v+, v- in (3. lOa, b), which measures the violation of equilibrium equations (3.2) must be

We consider the following decomposition of nodal displacements u = [I,, II,], where 8, are the degrees of freedom which represent the rigid body ‘displacements’ in the discretized model, and are taken to lie outside of the possible contact interfaces. Accordingly the following decomposition of (3.16a, b) is considered (A = [0, A])

604

G.

E.

STAVROULAKIS et al.

[z::::][::1+[A]=[L] (3.18a)

a2

WTZ = 0,

0,

(4. I b)

where

Au, + y = b.

(3.18b) * = [r, v+, v-IT,

Elimination of free variables II, between (3.18a) and (3.18b) result in the following LCP Rl&+Cr=Ik

(3.19a)

- Fr + GTdl + y = 6, ya 0,

r 20,

(3.19b) yTr =0,

(3.19c)

where R=K,,-K,&‘Kr,, C; = -K

P=P,-K,&‘PL.

K-IA~ ) I222

and 6, R are the appropriately transformed b and M (cf. [33]). We also note that the aforementioned interchange of the variables r and y does not affect the set of solutions of the LCP[34]. In order to solve LCP (4.1) we begin from a feasible starting point. This is not easy to obtain in problems with complex geometry. In order to overcome this difficulty, the algorithm which we use introduces a ‘measure of infeasibility’ expressed by an additional nonnegative variable z,. Problem (4.1) is thus written in the following form

F = AK;~IA~,

ti - ag - e,z, = 6,

6=b-ikKG’p2.

v? 2 0,

LCP (3.19) coincides with LCP (3.15), which was directly formulated previously based on a force method, if the free structure is isostatically supported by assuming ‘a, = 0’. Note that this formulation is more general than the force-like method due to the arbitrariness of the choice of 8, (and of the corresponding nodes). Concerning the interpretation of the discrete rigid body ‘displacements’ 1 used in (3.15) [respectively I, in (3.19)] we would like to mention that they actually represent the absolute displacements of the corresponding nodes. Obviously the values of the displacement of all nodes are dependent on the choice of the nodes 1 in (3.15) or of fi, in (3.19). 4. LEMKE’S ALGORITHM AND l-IS MECHANICAL INTERPRETATION

Matrix M of LCP (3.15) is positive semidefinite, due to the incorporation of the rigid displacements S in the formulation of the problem. In the general setting of the problem introduced in the previous section an equilibrium configuration may or may not exist, depending on the geometry of the whole structure (supported and unsupported parts) and the direction of the applied loading. If a solution does not exist equilibrium equations (3.2) are incompatible with the nonnegativity restriction of the unilateral contact forces (3.8a). Otherwise stated the set of admissible unilateral reaction forces r is empty. In order to avoid premature stopping of the algorithm, due to the presence of zeros in the diagonal of matrix M, we formulate the LCP (3.15) as follows (cf. analogous consideration in (32,331). By means of block principal pivoting on the principal submatrix F, + Fa of M we interchange the variable y with the variable r. LCP (3.15) is now written in the form *=KlRz+6,

2=[y,8+,B-lT

(4. la)

i 20,

(4.2a)

.ql> 0, WTZ= 0,

(4.2b)

where e,= [I, 1,. . . , IIT. We now are ready to describe the complementary pivot algorithm [35]. Let m be the number of variables contained in each vector 8, i in (4.2). Moreover, recall that an m-vector composed of certain nonzero elements of the vectors C and i is called a basic vector for the corresponding problem. Let vector * be the initial base vector of problem (4.2), thus i = 0. This means that initially all unilateral joints are closed (y = 0, respectively r > 0) and equilibrium equations are violated (v+ > 0, v- > 0). Let A be the m x (2m + 1) matrix and II the m-vector which are used to write equations (4.2a) in tableau form. Initially A and h take the forms [cf. (4.2a)] A = [II-R]

-e,,],

h = 6.

(4.3)

The algorithm has the following steps. Step I. Initialization. Set k = I. Let 8 be the initial basic vector. Find t” such that 6, = {min,h: j = 1, . , m}. If 6,, < 0 then continue with Step 2, otherwise the algorithm terminates with i: = 0 (unique if & > 0). Step 2. Initial infeasibility z&u.Drop from the basic vector the basic variable which corresponds to index rk and insert z. in its position. This is performed by going directly to Step 4. Step 3. Complementary change of basic variables. Set k = k + 1. Let s be the column of the nonbasic variable complementary to the dropped variable which has to enter the basis according to the complementary pivoting rule. Find the basic variable, to be

Rigid body displacements in unilateral problems

dropped, by means of the minimum ratio law, i.e. find tk such that R,, = min 3 : a,, > 0, j = 1, . . , m J

605

displacements or rotations. Then we get for the force method the LCP [cf. (3.15)] y-Fr-qz,=b,

,

U/J

y>

wherea,,,h,,i=l,..., 2m+l,j=l,..., m,arethe components of the current tableau A, h. If {a,,} < 0 for all j=l,... ,m, i.e. the algorithm is unable to introduce into the base the variable s without violating the nonnegativity restriction, then go to Step 6, otherwise continue with Step 4. This is done in order to drop from the basic vector the variable, which corresponds to index rk and insert in its position the nonbasic variable which corresponds to column s. Step 4. Apply the Gauss-Jordan pivot elimination step (first phase of the Simplex method of linear programming) with pivot row the r’th row and pivot column the sth column (respectively the (2m + I)th column] if k > I (respectively if k = 1). Step 5. Stopping criterion. If z,, is the dropped basic variable (i.e. z,, = 0) then go to Step 7, otherwise continue with Step 3. Step 6. Ray termination. The algorithm is terminated with ray termination (z,, > 0) and the LCP has no solution. Step 7. Normal exit. The algorithm is terminated with z0 = 0 and the current basic vector contains the solution of the LCP. From all the above steps we conclude that if ray termination occurs, i.e. z,, > 0, the LCP has no solution. In this case the data of the problem are incompatible, which means that the structure cannot be supported from the existing unilateral joints and for the applied loading and the geometry of the problem assumed. In particular, if ray termination occures one has that the set [cf. (3.2) and (3.8a)]

0,

r 2 0,

yTr = 0,

z,>O.

(4.7)

r-K’y-e,z,=b’, y 20,

r 20,

yTr= 0,

The infeasibility variable z,, expresses for the force method (respectively the displacement method) an additional interface opening (respectively prestressing force) which is applied to the unilateral contact interface in order to achieve displacement compatibility (respectively force equilibrium) in each step of the algorithm. We discuss first the force method (see Fig. 2a-c). The first value of z,, (i.e. zh’)) corresponds to an artificial opening required to separate all the possible contact nodes, except one pair of nodes which has an opening equal to zero. The subsequent steps of the algorithm gradully reduce the artificial opening z, from zk’)to zero by considering one-by-one the pairs of nodes which come into contact. Analogously in the displacement method (see Fig. 2a’-c’), zb’i corresponds to the value of a fictive prestressing force which is required in order to bring into contact all the possible pairs of contact nodes (see Fig. 2b’). The algorithm reduces the prestressing z0 from its initial value z&‘)to zero by deleting one-by-one the nodes where contact is not realized.

-s P

(4.4)

is infeasible (cf. [35]), which by Farka’s theorem (see e.g. [35, p. 5131) implies that the set

(4.6)

and for the displacement method

Gr=Bp,, r 2 0,

zoZ 0,

~$0 (a 1

r=O

(bl

z, = 0 (Cl

-y=iVG 0,

(4.5)

is feasible. This result will be compared with the conditions for the existence of a solution of the LCP which will be derived in the subsequent section. Alternatively one can use principal pivoting for the solution of the LCP with double pivoting within each step where a zero pivot appears in the pivoting operation (cf. the Cottle-Dantzig algorithm [36]). It is interesting to discuss the mechanical meaning of Lemke’s algorithm. Let us assume for the sake of simplicity, that here the bodies do not have rigid body

y=o

(a’)

( b’l

r&O

z, = 0 (C’l

Fig. 2. On the statical meaning of Lemke’s algorithm for unilateral contact problems.

G. E.

606

STAVROULAKU

et al.

5. EXISTENCE THEOREMS FOR UNILATERAL CONTACT PROBLEMS WITH PO!3lTVE SEMIDEFINITE MATRICES

We consider the unilateral contact problem [cf. the QPP (1.1) equivalently the variational inequality (l.S)]. In this section we shall give a condition which guarantees that the aforementioned problems with positive semidefinite matrix K have a solution. This condition is a special case of more general theorems proved for continuous variational inequalities by Fichera (see e.g. [1, 9). Let us consider first the QPP (1.1) [equivalently the variational inequality (1.5) formulation of the problem]. Let 9 be the space of rigid body displacements W= {UE R"lKu=O}.

(a)

v

(5.1)

Moreover let 51, = 9 n u,,

A

and32,={p/pEW,~-pE~,j. (5.2)

The following proposition gives a sufficient condition for the existence of a solution of problem (1.1) or (1.5) (with the general constraint Au - b < 0). Proposition 5.1. Suppose that either (i) 9, is bounded, or (ii) p’p < 0. for each p E W, -G&,

(b) Fig. 3. On the existence theorem.

Then by setting in (5.6) u* = u, + Lp, U,,E U,,, u*EUadasl++co,I,>O,Kp=O,pTp=O,wehave that

(5.3)

(Ku - P)~(u, + Ap - II) z o.

with Equivalently sup r:rsO,p’p#O,pEa,~E~, 1

= +a*

IIP II

(5.4)

and pTp = 0,

for each p E &,

(5.5)

then problem (1.1) or (1.5) has a solution u. Proof. The proof is a special case of a theorem of Fichera [I] where the constraint has the general form u E Dbwith K a convex closed set of a Hilbert space. In Fig. 3 we give an explication of the condition (5.4). Suppose that we consider on the plane the equilibrium of a material point submitted to x constraints Au - b ,< 0 which define the admissible set of Fig. 3a (BA is parallel to FA). Accordingly a material point at A (respectively at B or r) is in equilibrium, if it is subjected to a force P which forms an angle cp < 90” with A A (respectively BA or r A). The case of Fig. 3b is easier, since all directions p are of the type AA. (cf. in this context [1,371). Proposition 5.2. Suppose that (1.1) has a solution. Then for any p such that p E W, and satisfies (5.4) the condition pTp e 0 must be satisfied. Proof. If (1.1) has a solution u then it satisfies the variational inequality (1.5), i.e. CKu- P)T(n* - u) 2 0,

vu* E u,,,

(Ku- P)~(UO- n) _ pTp > o

I

u E u,.

(5.6)

,,

i.

which imples for A--, + co the condition (5.3) since p’p #O from (5.4). Q.E.D. Proposition 5.3. Suppose that II, is a solution of problem (1.1). Then any other solution u2 of the problem is given in the form u2 = UI + p,

where PEB,

p’p=O,

u,+pEu,d.

(5.7)

Both solutions give the same minimum. Proof. Suppose that u, E U,,,, u? E U,, are two solutions of (1.5). By setting u = u, , II* = u2 and u = u2, u+ = u, in (1.5), which is equivalent to (1.1) and adding we get

-04 -~2)~K(u,--u2)>0.

(5.8)

Inequality (5.8) implies due to the fact that K is positive semidefinite that p = u, - P*E 92. Then from and thus p’p ~0. (1.5) we get that p'(u,-u2)p0 Q.E.D.

Rigid body displacements in unilateral problems From the equivalent LCP formulation (3.15) the following proposition holds. Proposition 5.4. Let Fp > 0 for every p 2 0 such that Mp = 0. Then LCP (3.15) has a solution. Proof. The proof of proposition 5.4 results from proposition 5. I. Note in this respect that case (i) of proposition (5.1) is not possible for the LCP (3.15) and in case (ii) Rz = {O}. Q.E.D. Remark. The condition given in proposition 5.4. implies [cf. (3.1 SC)] that ~~6 = rTb - BTBp, > 0, for each p =[r,S+,S-]T>O such that Mp = 0 which implies that the set GTB < 0,

BTBp, < 0

is feasible. Thus the aforementioned condition coincides with the fact that the set (4.5) is infeasible, i.e. the set (4.4) is feasible and Lemke’s algorithm will never stop with ray termination. Existence theorems analogous to the one given previously have also been derived for general continuous variational or hemivariational inequality problems (in this respect see [2, 3, 38,391 and the references given there).

L

-35.0

-1 Fig. 4. First example, with one vertical rigid body displacement.

G. E.

608 6. NUMERICAL

STAVROULAKIS et al.

gaps arising at the frictionless contact interface as well as the values of the rigid body displacement calculated by means of the complementarity pivoting algorithm described before is given in Table 1 for the following material constants:

APPLICATIONS

The following applications reveal the advantages of the proposed method, i.e. the calculation of the rigid body displacements and rotations. (i) The unilateral contact problem between the parts of the structure shown in Fig. 4 is solved (one-half of the structure due to symmetry is shown and solved). The upper part of the structure has one rigid body displacement mode, the vertical displaceforce equilibrium ment 8,. The corresponding equation of the free part of the structure is written with respect to node 0 which is the pole mentioned in Section 3. A total number of 74 nodes and 65 finite elements are used for the finite element discretization of the structure. Triangular constant strain finite elements and rectangular isoparametric finite elements have been used as is shown in the figure. For a nodal vertical loading of magnitude p = 980 on node 0 the nodal contact forces and the respective

Case 1. Elastic moduli: EA = Es = 9.31 x 10’ (cf. the

elastic stamp problem). Case 2. Part A rigid, E, = 9.31 x 10’ (cf. the rigid

stamp problem). Case 3. EA = E, = 9.31 x 10’. Case 4. Part A rigid, EB = 9.31 x 10’.

Poisson’s ratio is taken to be v = 0.30 and the thickness of the plate is t = 1.0 (plane stress problem). In Fig. 5 the contact forces and the respective normal gaps at the interface are schematically depicted for cases 1 and 2. (ii) The second structure analysed is shown in Fig. 6. In this case the two displacement degrees of freedom of node 0 and the vertical displacement of + CASE 1 -=. CASE 2

1

2

3

4

6

6

7

6

9

10

11

INTERFACE NODES + CASE 1 0 CASE 2 1.2OE-01

T(b)

1 .OOE-Ol

6.OOE-02

:

6.00E-02

: 4.00E-02

2.OOE-02

O.OOE+OO 1

2

3

4

6

6

7

6

9

10

11

INTERFACE NODES

Fig. 5. Schematic representation of contact forces (a), and respective gaps (b), for the first example (cf. Table 1).

Rigid body displacements in unilateral problems

2x100.0

Y

3I

/ z’?

_

0

!

!

01

;

;

609

=f= $ <’ ‘!‘H

I

888888888888888 1+++++++++++++1 wwwwwwwwwwwwwww

G8888888888888Z ~ddddddddddddd~

L

35.0 -- -

-. ---_1

Fig. 6. Second example, with one vertical rigid body displacement and one rigid body rotation.

node 0, are used to simulate the three rigid body displacements’ of the free part A. For the finite element analysis 86 nodes and 77 finite elements are used. The number of the interface quantities calculated are given in Table 2, for the loading shown in Fig 6 and for the following material constants: v = 0.30, t = I.0 and the following combinations of elastic moduli: Case Case Case Case

I. E,=E,=9.31 x IO’. 2. Part A rigid, EB= 9.31 x 10’. 3. E,=E,=9.31 x IO’. 4.Part A rigid, E, = 9.31 x IO’.

The results of cases 1 and 2 are given schematically in Fig. 7. (iii) For comparison we also solve the same example by adding two additional soft springs at the nodes Z and H (Fig. 6) in order to eliminate the rigid body displacement modes. In this example we have taken the stiffness of each additional spring to be equal to 1.61% of the value of the corresponding (diagonal) member in the global stiffness matrix of the structure. This low additional stiffness is enough to eliminate the rigid body ‘displacements’ of the structure which can be solved in the sequel by means of more simple algorithms of unilateral contact analysis. The results given in Table 3 are obtained by means of an appropriate application of the Hildreth and d’Esopo iterative method (cf. [28]) for the analysis of the QPP (1.5) with positive definite

668888888888866 Ill+++++++++lll wwwwwwwwwwwwwww

~~=888888888=~~ tiA_id&&d&&&&&+Ati 8888 Z~SZ~ZZ8888 +++++++++++++++ wwwwwwwwwwwwwww

888~z?~%~%Zz36888 &~~.&Apjpjgl+&Auj~~~ %8888888888888?? 1+++++++++++++1 wwwwwwwwwwwwwww 28888888888888% &dddddddddddddG

82 II ww

“F A+

610

G.

E.

STAVROULAKIS et al. -CASE1

3mE+91

=-CASE2

T(a)

3.WE+91 : N 2AOE+Ol 1 t ?..WE+Ol T F l.!iOE+Ol 0 ! :

l.WEtOl 5.WE+O9

o.ooE+oo 1

2

3

4

5

6

9

7

9

10

11

12

13

14

15

INTERFACE NODES *CASE1

1

2

3

4

5

6

7

eCASE2

9

9

10

11

12

13

14

15

INTERFACE NODES

Fig. 7. Schematic representation of contact forces (a) and respective gaps (b) for the second example (cf. Table 2).

Table 3. Second example, case I: comparison between the contact forces obtained through the fictive soft sorinn conceot and the exact method DroDosed in this oawr Nodal contact forces Difference Direct method D=C-B C

Interface pair of nodes A

Fictive spring B

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

O.OOE+ 00 5.4OE - 01 6.30E + 00 1.20E + 01 1.70E+Ol 2.20E + 01 2.80E + 01 2.70E + 01 2.80E + 01 2.20E + 01 1.70E + 01 1.20E+Ol 6.4OE + 00 5.4OE - 01 O.OOE+ 00

O.OOE+ 00 2.288 + 00 7.30E + 00 1.24E+Ol 1.72E+Ol 2.13E+Ol 2.68E + 01 2.548+01 2.688 + 01 2.13E+Ol 1.72E + 01 1.24E+Ol 7.23E + 00 2.28E + 00 O.lWE+ 00

O.OOE+ 00 1.74E + 00 7.70E - 01 4.OOE- 01 2.OOE- 01 7.OOE- 01 1.20E + 00 16OE+OO 1.20E+OO 7.OOE- 01 2.OOE- 01 4.OOE- 01 7.70E - 01 1.74E + 00 O.OOE+ 00

76.31% 13.41% 3.22% 1.16% 3.28% 4.47% 6.29% 4.47% 3.28% 1.16% 3.22% 13.41% 76.31%

Total force

198.68

200

1.32

0.66%

Deviation E = (D/C) x 100

Rigid body displacements in unilateral problems +

c

3.00E+Ol

T

25OE+Ol

--

FICTIVE SPRING

611

-w DIRECT METHOD

0 N T 2.OOE+Ol -A

1

2

3

4

9

9

7

9

9

10

11

12

13

14

15

INTERFACE NODES

Fig. 8. Comparison of the contact forces of the second example calculated by means of the fictitious soft springs concept and by using the here proposed method (exact solution cf. Table 3).

Table 4. Interface results for the structure of Fig. 6 under

matrix due to the presence of the additional springs. Moreover contact forces of Table 3 are summed up to give a total vertical reaction equal to 198.68 while the applied vertical load is equal to 200. Obviously the force which is missing from equilibrium (equal to 1.32 or 0.66% of the external loading) corresponds to the forces of the fictive additional springs. Moreover the comparison of the contact forces distribution between exact and approximate solution (Table 3) shows the insignificant influence of these additional fictive springs. Analogous observations concern the stress and strain distributions in the structure. (iv) The structure of the second example (with EA = EB = 9.3 1 x 10’) is also solved with two downward loads of magnitude 100 placed on the nodes 2 and H (cf. Fig. 6). The results summarized in Table 4 and Fig. 9 show the different behaviour at the interface in this case.

vertical loading at the points Z and Interface pair of nodes

Contact force

I 2 3 4 5 6 7 8 9 10 II I2 I3 14 I5

2.27E + 01 3.28E + 01 2.20E + 01 1.47E + 01 5.14E+OO 2.59E + 00 OBOE+ 00 OBOE+ 00 OBOE+ 00 2.59E + 00 5.14E + 00 I .47E + 01 2.20E + 01 3.28E + 01 2.27E + 01

Displacement of pole 0

2.4lE - 05

Cap OBOE+ OBOE+ O.OOE+ O.OOE+ OBOE+ O.OOE+ 3.20E 2.66E 3.20E O.OOE+ O.OOE+ O.OOE+ OBOE+ OBOE+ O.OOE+

LOADING AT 2 AND H

3.9OE+Ol

;

2.00E+Ol

T F

l.mE+Ol

ii c

l.OOE+Ol

E &OOE+OO

1

2

3

4

5

6

7

9

9

10

11

12

13

INTERFACE NODES

Fig. 9. Schematic presentation of the contact forces given in Table 4.

14

H

lb

00 00 00 00 00 00 08 08 08 00 00 00 00 00 00

G. E. STAVROIJLAKISet al.

612

Fig. 10. Two-dimensional rigid object in a multifingered elastic

(v) We now give an application of the theory to the grasping of a rigid object by means of an elastic

multifingered robotic hand. Let us consider the twodimensional rigid body shown in Fig. 10, which is grasped by means of six elastic potential fingers with flexibilities Fi = 0.1, i = 1,. . . ,6, with zero initial fingertipobject gap. Let the whole system be referred to the orthogonal Cartesian coordinate system Oxy. The fingertip contact forces and the corresponding openings (gaps) as well as the rigid body displacements and rotations required in order to reach the object its equilibrium configuration are given in Tables 5 and 6 for the following sets of externally applied loading respectively: Table 5: loading equal to 10 exerted from finger 1; Table 6: Px = -30, Py = - 20, M = - 9. The fingertip contact forces and

the applied loading constitute what is called in robotics a ‘force closure’ system for the rigid object [40]. In order to investigate further the assumed gripper problem we apply a loading of the following form: Px = P cos fJ cos cp, Py = P cos 0 sin cp,M = P sin 0. Obviously by assuming all combinations between O,
Table 5. Solution of the elastic gripper shown in Fig. 10 under the external loading {Px = -30, Py = -20, M = -9} on the rigid object

Finger 1 2 3 4 5 6

Contact reaction and corresponding gap Contact Fingertip reaction displacement Gap 0 0.296D + 01 0 0 0.858D + 00 0 O.l16D+02 0 o.l16D+ol 0.274D + 02 0 0.274D + 01 0 0.162D +01 0 0.134D + 02 0 0.134DfOl

Rigid body displacements Displacement Displacement component value n: II+ Qd+ fi; 5 cp-

0 o.l61D+OO 0.949D - 01 0.365D + 01 0 0

Table 6. Solution of the elastic grinner problem shown in Fig. 10 under the assumption that a force of magnitude IO-is applied in the position of finger 1 _

Finger 1 2 3 4 5 6

Contact reaction and corresponding gap Fingertip Contact displacement reaction Gap 0 0.213D + 01 0 0.752D + 01 0 0.752D + 00 0.752D + 01 0 0.752D + 00 O.lOOD+02 0 O.lOOD+ 01 0 0.693D+OO 0 0 0.693D+OO 0

Rigid body displacements Displacement Displacement value component fi: V ‘p+ h,_ 4 Q-

0 0.8OOD+ 00 0.915D - 01 0.143D + 01 0 0

Rigid body displacements in unilateral problems

613

3

4

+=/2

n

h/2

2n

(a)

(b) Fig. 11. Complete investigation of the gripper-object system shown as shown in (b). REFERENCES

1. G. Fichera. Boundary value problems in elasticity with unilateral constraints. In Encyclopedia of Physics (Edited by S. Fliigge), Vol. VI, a/2. Springer, Berlin (1972). 2. G. Duvaut and J. L. Lions, Les Inhpnfions en Mknique et en Physique. Dunod, Paris (1972). 3. P. D. Panagiotopoulos, Inequality Problems in Mechanics and Applications. Convex and Nonconvex Energy Fundions. BirkhHuser, Base1 (1985). (Russian trans-

4. 5. 6.

7.

lation, MIR, Moscow, 1989.) I. Hlavafek. J. Haslinger, J. N&s and J. LoviSek, Solution of Variational Inequalities in Mechanics. Springer, New York (1988). P. D. Panagiotopoulos, Non-convex energy functions. Hemivariational inqualities and substationarity principles. Acfu Mechnnicu 48, 1I l-1 30 (1983). P. D. Panagiotopoulos, Hemivariational inequalities and their applications. In Topics in Nonsmoorh Mechunits (Edited by J. J. Moreau, P. D. Panagiotopoulos and G. Strang), pp. 77-142, Birkhauser. Base1 (1988). P. D. Panagiotopoulos, Nonconvex superpotentials and hemivariational inqualities. Quasidifferentiability in

8.

9.

10.

Il.

12.

in

Fig. 10. Active fingers are denoted

mechanics. In Nonsmooth Mechanics and Applications. CISM Leerwe Notes Vol. 302 (Edited by J. J. Moreau and P. D. Panagiotopoulos), pp. 83-176. Springer, Vienna (1988). P. D. Panagiotopoulos and G. E. Stavroulakis, A variational hemivariational inequality approach to the laminated plate theory under subdifferential boundary conditions. Q. Appl. Muth. XLVI, 409430 (1988). P. D. Panagiotopoulos. Semicoercive hemivariational inequalities. On the delamination of composite plates. Q. Appt. Math. XLVl, 61 l-629 (1989). P. D. Panagiotopoulos and D. Talaslidis, A linear analysis approach to the solution of certain classes of variational inequality problems in structural analysis. inr. J. Solidr Strucr. 16, 991-1005 (1980). D. Talaslidis and P. D. Panagiotopoulos, Linear finite element approach to unilateral contact problem in structural dynamics, Int. J. Numer. Mefh. Engng 18, 1505-1520 (1982). E. Mitsopoulou and I. N. Doudoumis, A contribution to the analysis of unilateral contact problems with friction. Solid Mech. Arch. 12, 165-186 (1987).

G. E. STAVROULAIUS et al.

614

13. C. D. Bisbos, A Cholesky condensation method for unilateral contact problems. Solid Mech. Arch. 11, l-23 (1986). 14. G. Nitsiotas, Die Berechnung statisch unbestimmter Tranwerke mit einseitiaen Bindunnen. Archiv _ Innenieur _ 41, &-6o (1970). 15. G. Maier, Mathematical programming methods in structural analysis. In Variational Methods in Engineering (Edited by C. Brebbia et al.), Vol. II, 8/l. Southampton University Press, Southampton (1973). 16. J. J. Moreau, La notion de sur-potentiel et les liaisons unilaterales en Clastostatique. C. R. Acad. Sci. Paris

26. J. L. Lions, Optimal Control of Systems Governed by Partial Differential Equafions. Springer, Berlin (1971). 27. Z. P. BaZant, Imbricate continuum and its variational derivation. J. Engng Mech. ASCE 110, 1693-1712 (1984). 28. H. P. Kiinzi and W. Krelle, Nonlinear Programming. Blaisdell (1966). (German edition, Springer,l962.) 29. M. M. Kostreva. Generalization of Murtv’s direct algorithm to linear and convex quadratic programming.

267A, 954-957 (1968). 17. C. D. Bisbos, A discrete nonlinear boundary variational

Ingenieur Archiv 52, 63-76 (1982). 3 1. H. Bufler, Derivation of the variational inequalities and

inequality for 3-D frictional contact problems. Internal Report 1542, Institute of Steel Structures, Aristotle University, Thessaloniki (1990). 18. P. D. Panagiotopoulos, A nonlinear programming approach to the unilateral contact and friction boundary value problem in the theory of elasticity. Ingenieur

extremun principles of the frictionless elastic contact problem. Comp. Meth. Appl. Mech. Engng 53, 163-182 (1985). 32. L. Johanson and A. Klarbring, The rigid punch problem with friction using variational inequalities and linear complementarity. In L. Johanson, Rigid punch problems and Green’s functions, Thesis No. 204, Dept of Mechanical Engineering, Linkiiping Institute of Technology, Sweden (1990). 33. A. A. Cannarozzi, On the resolution of some unilaterally constrained problems in structural analysis. Comp.

Archiv 44, 421432

(1975).

19. I. N. Doudoumis and E. N. Mitsopoulou, On the solution of the unilateral contact frictional problem for general static loading conditions. Cornput. Struct. 30, 111l-l 126 (1988). 20. M. Z. Cohn and G. Maier (Eds), Engineering plasticity by mathematical programming. Proceedings of the NATO Advanced Study Institute, University of Waterloo 1977. Pergamon Press, New York (1979). 21. P. D. Panagiotopoulos, Stress-unilateral analysis of discretized cable and membrane structures in the presence of large displacements. Ingenieur Archiv 44, 291-300 (1975).

22. 0. Z. Zienkiewicz, S. Valliapan and I. P. King, Stress analysis of rock as a ‘notension’ material. Geotechnique 18, 56-66 (1968).

23. G. Maier, A quadratic programming approach to certain classes of nonlinear structural problems. Meccanica 3, 121-130 (1968).

24. M. Fremond, Etude de structures viscoelastiques stratifiees soumises a des charges harmoniques et de solides Clastiques reposant sur ces structures. These de doctorat d’etat, Univ. Paris VI (1971). 25. I. Ekeland and R. Temam, Convex Analysis and Variational Problems. North-Holland, Amsterdam (1976).

J. Optimiz. Theory Applic. 62, 63-76 (1989). 30. H. Bufler, H. Lieb and G. Meier, Frictionless contact

between an elastic stamp and an elastic foundation.

Meth. Appl. Mech. Engng 24, 339-357 (1980). 34. K. G. Murty, On the number of solutions

to the complementarity problem and spanning properties of complementary cones. Linear Algebra Applic. 5,65-108

(1972). 35. K. G. Murty, Linear Complementarity. Linear and Nonlinear Programming. Heldermann, Berlin (1988).

36. G. B. Dantzig and R. W. Cottle, Positive semi-definite programming. ORC 63-18 (RR), Operations Research Center, University of California, Berkeley (1963). 37. J. L. Lions and G. Stampacchia, Variational inequalities. Commun. Pure Appl. Math. XX, 493-519 (1967). 38. J. J. Moreau, P. D. Panagiotopoulos and G. Strang (Eds), Topics in Nonsmooth Mechanics. Birkhiiuser, Base1 (1988). 39. J. J. Moreau and P. D. Panagiotopoulos (Eds), Nonsmooth Mechanics and Applications. CISM Notes Vol. 302. Springer, Vienna (1988).

Lecture

40. J. Kerr and B. Roth, Analysis of multifingered hands. Int. J. Robotics Res. 4, 3-17 (1986).